Generating entangled quantum microwaves in a Josephson-photonics device

When connecting a voltage-biased Josephson junction in series to several microwave cavities, a Cooper-pair current across the junction gives rise to a continuous emission of strongly correlated photons into the cavity modes. Tuning the bias voltage to the resonance where a single Cooper pair provides the energy to create an additional photon in each of the cavities, we demonstrate the entangling nature of these creation processes by simple witnesses in terms of experimentally accessible observables. To characterize the entanglement properties of the such created quantum states of light to the fullest possible extent, we then proceed to more elaborate entanglement criteria based on the knowledge of the full density matrix and provide a detailed study of bi- and multipartite entanglement. In particular, we illustrate how due to the relatively simple design of these circuits changes of experimental parameters allow one to access a wide variety of entangled states differing, e.g., in the number of entangled parties or the dimension of state space. Such devices, besides their promising potential to act as a highly versatile source of entangled quantum microwaves, may thus represent an excellent natural testbed for classification and quantification schemes developed in quantum information theory.


I. INTRODUCTION
The concept of entanglement is at the heart of quantum physics: as one of the cornerstones at the foundation of quantum mechanics and, at the same time, as a key ingredient in emerging quantum technologies introducing a new quantum resource into communication, computing, and sensing.
These two interconnected aspects of entanglement are reflected in a branching of research interest. On one hand, a more abstract, mathematical direction of quantum information theory which maps out the boundary between classical and quantum world and further charts quantum states into various classes of entanglement. More and more sophisticated (qualitative) witnesses and (quantitative) measures of entanglement are devised in this field [1][2][3][4][5][6][7][8] with many open questions remaining, in particular, concerning mixed states and multipartite systems [9,10]. Oftentimes research will consider complex quantum states characterized by certain, special classes of density matrices, for which entanglement witnesses or measures can be worked out, while typically the question in which actual systems and under which actual circumstances such states may naturally be encountered garners less attention.
On the other hand, proposed and realized entanglement sources in the optical and microwave regime typically strive to realize particularly simple entangled states, for instance, pure (squeezed) Gaussian states with small, negligible admixtures [11][12][13], NOON states [14][15][16], or maximally entangled states [17][18][19], thereby accessing only a small number of all interesting types of entangled states. In the microwave regime, besides setups using analogues of typical nonlinear optical elements [13], more recent schemes extend the ability to create arbitrary quantum states in a microwave resonator [20] to several, spatially separated modes [15]. While by properly choosing an elaborate pulse scheme, in principle, any entangled state could be created, it is typically a maximally entangled or another simple pure entangled state which such experiments aim for.
The creation of multipartite or other more complex entangled states requires increasingly complex pulse schemes. These are reachable in such systems since one can build on the immense research effort (and the resulting amazing progress in performance and control) which has been spent on these standard circuit quantum electrodynamic (QED) setups as part of the larger quest for universal quantum computing. Here we argue, however, that combining two key elements of circuit-QED setups, namely, Josephson junctions and microwave cavities, in a much simpler, less demanding device can offer an alternative, fully tunable and versatile entanglement-generating source.
The nontrivial entangled states, which appear as (quasi-)stationary states of these voltage-driven, damped systems, allow to explore the richness of entanglement phenomena along several dimensions: the number of entangled parties, continuous or discrete degrees of freedom, fully entangled or biseparable states, all this is controllable in a single experimental device. In turn, the full wealth of entanglement witnesses and measures has to be used in attempt to characterize the entanglement properties of these versatile entanglement sources to the fullest possible extent.
Josephson-photonics devices or other entanglement  sources utilizing similar entanglement creation mechanisms, hence, highlight the necessity of better bridging the gap between the abstract quantum-information branch of research on entanglement and more practicaloriented approaches aimed at utilizing entanglement as a resource for various quantum technological applications.
In this paper, we will first briefly discuss a basic example of realizing either bi-or tripartite entanglement by a simple change of the voltage bias applied to the Josephson junction (Sec. II). After having thus demonstrated that such a device can, in principle, be used as a source of genuinely multipartite entangled microwaves, we will explain how also other characteristics of entanglement can be explored by varying other experimental "knobs". A detailed discussion of bi-and multipartite entanglement in Secs. III and IV, respectively, then reveals the full richness of entanglement phenomena in Josephsonphotonics devices. Limitations which occur in current experimental setups and possible improvements and modifications which might thus become necessary are discussed in Sec. V. We conclude in Sec. VI and address open questions which remain for future research.
Biasing the junction so that its Josephson frequency ω J = 2eV / matches the sum of the modes frequencies, ω J = q ω q , the transfer of a single Cooper pair across the Josephson junction gives rise to a simultaneous creation (or annihilation) of one photon within each of the modes. This common creation process is balanced by subsequent individual leakage of photons via output lines connected to each of the different resonators so that eventually a stationary [42] (and possibly entangled) state of the spatially separated cavities is reached.
Formally, the particular resonant processes at the given bias are picked out by a rotating-wave approximation after the system Hamiltonian is transformed to a frame rotating with the Josephson frequency ω J [28,29,32,33]. The resulting effective time-independent Hamiltonian contains a renormalized Josephson energy, E * J = E J q √ κ q e −κq/2 , and normal-ordered Bessel functions J 1 of the first kind. These are functions of the photonic number operators n q = a † q a q of the various resonators, expressed in terms of creation/annihilation operators with [a q , a † q ] = 1, and the dimensionless parameter κ q = (2e 2 / ) L q /C q measuring each oscillator's zeropoint quantum fluctuations. In the Appendix A, we discuss the equivalent structure of the Hamiltonian for other resonances, e.g., if the chosen bias allows single photon excitations in some subset of all coupled oscillators.
The bracketed term in Eq. (1) indicates the fundamental photon creation/absorption process considered here: each forward (backward) tunneling event of a Cooper pair goes along with the simultaneous creation q a † q (absorption q a q ) of one photon in each of the oscillator modes. These fundamental processes, however, are modified by the inherent nonlinearity of the Josephson junction, reflected in the appearance of the product of Bessel functions, which makes the driving nonlinear.
Generally speaking, the effect of these nonlinearities on the dynamics of the system is governed by two different experimental parameters, the κ q parameter(s) and the Josephson energy E J . The former fixes the relative size of the transition matrix elements between neighboring states of oscillator q and thus effectively sets its level structure as we will see below in Sec. III B. The latter determines the driving strength and thus the population of the corresponding levels. As we will see later on [cf. Fig. 1(c)], varying these parameters (and the resonance chosen) allows to explore a wide range of entanglement phenomena in Josephson-photonics devices.
Creating and annihilating photons according to the Hamiltonian in Eq. (1), a stationary state is reached in due course as excited photons leak out of the resonators into the electromagnetic environment after a lifetime 1/γ q . Focusing on the zero-temperature limit, this dynamics can be described by a quantum master equation of the standard Lindblad form [43] for the density operator ρ of the cavity degrees of freedom.
A first basic example demonstrating the entanglement power of the device and at the same time illustrating the general strategy pursued in our investigation is sketched in Fig. 1(b). Let us consider a setup consisting of three cavities a, b, and c, with the photonic operators now denoted a † 1 = a † , a † 2 = b † , and a † 3 = c † . Selecting two biasing conditions ω J = 2eV / = ω a + ω b (+ω c ) and consequently two different resonant excitation processes, the driving is turned up to investigate if (and what type of) entanglement is achieved.
What statements on entanglement can be made depends on the choice of entanglement witness or measure. Later on, we will present a number of more elaborate witnesses designed to optimize the information on entanglement gained. For the moment, however, we consider a particularly simple witness experimentally accessible via measurements on the output lines. Namely, the violation of the inequalities [44,45] | a † b † c † st | ≤ n a n b st n c st + perm.

sym.
= 3 n a n b st n c st (4) bears witness to the presence of genuine bi-and tripartite entanglement for the full range of driving shown in Fig. 1(b). The equality signs hold for the symmetric case, κ a = κ b (= κ c ) and γ a = γ b (= γ c ), also assumed in the figure.
Indeed, we show data for small κ a/b/c = 0.1 so that nonlinearities appear for large occupation numbers only. In the weak driving case, there are but small corrections to the bare entangling process of a nondegenerate parametric amplifier with driving ∝ E J a † b † + c. c and its three-party equivalent. In consequence, for the ω J = ω a + ω b resonance we immediately recognize in Fig. 1(b) the known parametric-amplifier results, n a st = n b st ∝ E 2 J and a † b † st ∝ E J for weak driving [32,35]. If the inequality holds, no statement on entanglement is possible based on the chosen witness. Similarly, for the ω J = ω a + ω b + ω c resonance we cannot witness entanglement between two modes alone ( a † b † st ≡ 0), but instead all three resonators are entangled for the whole driving range shown.
These simple cases can now be used to exemplify the key strategy of this investigation. The versatile entangling properties of Josephson-photonics devices are exploited to explore the space of entanglement phenomena along several directions. Schematically, this is visualized in Fig. 1(c). Roughly speaking, different experimental knobs in our setup correspond to different dimensions of the complete space of all possible entangled states: (i) Changing the bias voltage, the number of actively involved oscillators and thereby the number of parties which are (potentially) entangled is varied. (ii) Changing the parameter κ q , different effective Hilbert (sub-)spaces for party q are realized. For κ q → 0, for instance, all levels of oscillator q are accessible, while for special values of κ q oscillator q reduces to an effective Nlevel state (most importantly κ q = 2 yields a two-level system as discussed in Sec. III B below). (iii) The driving strength E J finally is changed to tune the system through different classes of entanglement. These can be as simple as in the basic example above, where our witness only allows the identification of an entangled region below some E J and an unchartered region above. The aim, however, (achieved in this paper for the cases marked by the blue and red regions in the figure) is the complete characterization of entanglement for the chosen scenario.

III. BIPARTITE ENTANGLEMENT
In the three-oscillator setup presented above, a bipartite system is realized by selecting a bias condition 2eV / = ω a + ω b . The corresponding effective Hamiltonian differs from an experimentally already realized twocavity setup merely in a renormalization of the Josephson energy, see Eq. (A1) in the Appendix A. In any case, the simultaneous creation of a single photon in each of the oscillators a and b leads to strongly correlated dynamics  including potentially, but not necessarily, quantum entanglement. These two subsystems are called entangled if their corresponding mixed state ρ is not separable, i.e., if it can not be written as a convex combination of product states of the subsystems: ρ = j p j ρ j a ⊗ ρ j b with convex weights p j > 0 and j p j = 1.
The question whether a given bipartite density matrix is separable or entangled is still lacking a general answer. Over the last decades, however, a large number of different criteria have been proposed to detect the presence of entanglement or to prove separability (for an overview see, e.g., Refs. [1][2][3][4][5][6]). In experiments, the analysis of entanglement is often restricted to simple witnesses based on directly measurable quantities [see, e.g., Eq. (3) in the previous section]. Our theoretical approach, however, provides full information on the system, i.e., the complete density matrix, which is experimentally accessible only by full state tomography, so that in essence the whole range of theoretically established entanglement criteria may be applied.

A. PPT criterion and logarithmic negativity
An entanglement criterion both powerful and simple is the criterion based on positive partial transposition (PPT criterion) [46]. It relies on the fact that taking the transposition with respect to a single subsystem only does not necessarily map a state ρ onto a quantum state again. If, however, ρ is a separable state, its partial transpose ρ Ta = j (ρ j a ) T ⊗ ρ j b (and analogous for subsystem b) indeed represents a valid density matrix and is thus positive semidefinite, ρ Ta ≥ 0 ⇔ ρ T b ≥ 0. In consequence, we can conclude that ρ is entangled if its partial transpose is not positive semidefinite, i.e., if it has at least one negative eigenvalue. In this case, we say that ρ has a negative partial transpose (NPT). Importantly, the converse is explicitly not true in general: a PPT, i.e., the partial transpose of ρ is positive semidefinite, does only imply separability for the special cases of the low-dimensional 3 × 2 and 2 × 2 systems [47,48].
To additionally quantify the amount of entanglement of ρ, the logarithmic negativity E N [ρ] = log 2 (||ρ Ta || 1 ) [49,50], which is directly linked to violation of the PPT criterion, is used. This entanglement measure is based on the trace norm || . . . || 1 of the partial transpose ρ Ta , which is related to the sum of the negative eigenvalues λ j − of ρ Ta : ||ρ Ta || 1 = 1 + 2| j λ j − |. PPT criterion as well as logarithmic negativity are easy to calculate for arbitrary-dimensional state spaces if the density matrix is known; however, they obviously fail in detecting some entangled states, namely, those with a positive partial transpose.
We can now proceed to pursue our general strategy and explore the remaining directions of the space of entanglement phenomena (sketched in Fig. 1(c) as a plane) for this two-party case of Josephson photonics. For a first overview, the symmetric κ q direction, i.e., κ a = κ b = κ, is chosen and also the damping of the two oscillators is assumed to be equal, γ a = γ b = γ. Fig. 2(a) then shows the driving dependence of logarithmic negativity E N , while κ is allowed to increase from the parametric-amplifier limit κ → 0 with its harmonic level structure and driving. The results for the entanglement witness shown in Fig. 1(b) are close to this limit with κ = 0.1.
Two important features of Fig. 2(a) will be explored in the following discussion. Firstly, the logarithmic negativity E N shows a maximum which becomes the more pronounced the smaller the value of κ is with its position monotonically shifting from At κ = 2, however, E N does not only take a minimum but vanishes and stays zero for all 2 ≤ E J /E c2 J . This feature, specific to κ = 2, is related to the restricted Hilbert space of each of the entangled parties, as the suppressions of any transitions to higher excited levels effectively reduces the harmonic oscillator to a two-level system.

B. Restricted Hilbert space
In fact, the dashed lines in Fig. 2(a) [corresponding to the cross sections in Fig. 2(b)] indicate those special values of κ where the state space of each of the two oscillators is effectively reduced to a sixteen-(κ ≈ 0.23), three-(κ = 3 − √ 3), and two-level system (κ = 2). These special values simply follow from the roots of the transition matrix elements T ma, where the nonvanishing matrix elements are a function of E J and γ.
5) (the inverse of the golden ratio). The entanglement features of the 2 × 2 system are passed on nearly unchanged to systems with unequal κ q chosen such that one party still is a two-level system, while the dimension of the Hilbert space of the second party is increased, e.g., the 16×2 realization in Fig. 2(b). The similarity to the 2 × 2 case can be explained by the fact that, in consequence of energy conservation and the simultaneous creation of photons, the mean occupations of the two oscillators are coupled, γ a n a st = γ b n b st . The presence of the qubit thus leads to signatures of state-space restriction in the sixteen-level system (cf. Refs. [35,39]), which are reduced (γ a < γ b ) or enhanced (γ a > γ b ) for asymmetric decay rates.
Note that E N [ρ st ] = 0 does, in general, not imply separability for a 16 × 2 system. However, we know that if an N × 2 state ρ is invariant under partial transposition with respect to the qubit, i.e., ρ T b = ρ, then ρ is separable [51]. Since the property of separability is independent of invertible local transformations U acting on C 2 , also off with ρ diag and ρ off being the diagonal and off-diagonal part of ρ st , respectively. Choosing Separability for strong driving in the N × 2 cases does not mean that quantum correlations between the two subsystems are absent. This is reflected in a nonzero quantum discord [52,53], which measures quantum correlations by the difference between quantum mutual information and classical correlations. States of the 2 × 2 and 16 × 2 system with E N [ρ st ] = 0 have a strictly positive quantum discord, i.e., they contain nonclassical correlations even though they are separable. This is shown using a simple witness [54] which allows the detection of nonvanishing quantum discord: If a bipartite state ρ does not commute with ρ a ⊗ 1, where ρ a = Tr b {ρ}, then its quantum discord is nonzero and ρ is nonclassically correlated.

C. Modified driving Hamiltonian
Steady-state entanglement vanishes for strong driving only in case that a qubit is involved. For all other realizations, e.g., the 3×3 or 16×16 realization, the logarithmic negativity takes a minimum around E J /E c2 J = 2 followed by an increase in E N . This increase can be traced back to coherences associated with higher photon creation processes, e.g., ρ 00,22 in the 3 × 3 case. These do not vanish but saturate in the limit E J /E c2 J → ∞. If the absolute value of these coherences in this limit is sufficiently large compared to the populations, entanglement can occur even in the regime of very strong driving.
In this section, we want to focus on this strong-driving limit to illustrate that the κ direction of our sketch of entanglement space above is more complex than suggested so far by merely discussing the reduced Hilbert space of the two parties. To that end, we restrict the cavities to a 3 × 3 system but study the impact of modifying the relative size of the different transition matrix elements of the driving Hamiltonian. This is not purely an academic exercise but is motivated by recent circuit-QED experiments which employ additional oscillating driving fields to essentially cut off the harmonic ladder of energy states of a stripline cavity at some level N [55]. We envision similar schemes being accomplished in our Josephson-photonics setup and thus study a so-realized 3 × 3 system while allowing for a variable κ to change the ratio T 1,2 /T 0,1 of the single-cavity matrix elements, see Fig. 3.
The dashed lines indicate here some special realizations of 3 × 3 systems: (i) an otherwise harmonic system (κ → 0) cut off to 3 × 3, (ii) a symmetric 3 × 3 system (κ = 2− √ 2), (iii) the "native" Josephson-photonics 3×3 system (κ = 3 − √ 3), and (iv) the 2 × 2 system (κ = 2) also natively realized in Josephson photonics without additional fields. The results depicted in Fig. 3 show that whether ρ st in the strong-driving limit is entangled or not strongly depends on the value of κ and thus on the ratio of the matrix elements. Among the indicated special cases, only the native Josephson-photonics 3 × 3 system is entangled in the strong-driving limit.

D. Entanglement dynamics
The time-resolved statistics of photon emission events from the cavities has recently been studied experimentally and theoretically for Josephson-photonics setups [25][26][27][33][34][35][36][37][38]. Borrowing tools from quantum optics [56][57][58], in particular the second-order correlation function g (2) (τ ), provides a deep insight into the strongly correlated quantum dynamics of the oscillator subsystems. Following this work, it seems quite natural to also consider the entanglement dynamics after a single photon emission event is observed. The photon emission process from oscillator a is described by a jump operator J a acting on the density operator, J a = γ a aρa † .
The dynamics of entanglement is also often studied in the context of decoherence. Then, an initially prepared entangled state interacts with an environment leading to a decay and (oftentimes) a sudden death of entanglement after a finite time [59][60][61][62].
Here, we consider again the special case of a 16 × 2 system moderately driven, E J /E c2 J = 1.75, below the entanglement-separability threshold to its steady state.
Then at τ = 0 − , a photon leaving the system is detected. Figure 4 pictures the subsequent time evolution of the logarithmic negativity after emission from cavity a or b for γ a /γ b = 1 (solid lines) and γ a /γ b = 0.3 (dashed lines). In all these cases, the two oscillator subsystems are separable at τ = 0 + immediately after the photon emission event. Note that E N = 0 here is sufficient for separability after a photon measurement, which can be proven by a similar reasoning as outlined above in Sec. III B for the steady state of the 16 × 2 system. If the emitted photon stems from the two-level system leaving it in its ground state, separability is obvious. Surprisingly, however, a separable state is also found if the photon is emitted from the sixteen-level system.
Let us first focus on symmetric decay rates, i.e., γ a /γ b = 1. If a photon jump out of the sixteen-level system is detected, the two-level system might still be in its excited state blocking further photon creation processes and thus hindering the development of entanglement. This is why E N stays zero for some time after the emission event during which it becomes more and more likely that the two-level system has already relaxed before bipartite entanglement is created again. In contrast, if a photon is emitted from the two-level system, a new creation process can occur immediately, which is reflected in a strong increase in E N on short time scales. In this case, the entanglement features a distinct maximum but vanishes again for a finite span of time before E N finally approaches its steady-state value. Note that the "suddenness" of death and revival features, here as well as in the decoherence dynamics of entanglement, is of course simply related to the definition of logarithmic negativity: while all density matrix elements and eigenvalues evolve smoothly, the sum of negative eigenvalues only does not.  Figure 3. Steady-state logarithmic negativity EN in the strong-driving limit, EJ /E c2 J → ∞. State space is cut to a 3 × 3 system of two symmetric oscillators. Whether steadystate entanglement is observable depends on the κ parameter, which determines the ratio T1,2/T0,1 of the two transition matrix elements. Dashed lines indicate special realizations of 3 × 3 systems (see main text for details).  Figure 4. Dynamics after a photon measurement of the logarithmic negativity EN for a 16 × 2 system initially in steady state for EJ /E c2 J = 1.75. Observing a photon leaving oscillator a (magenta) or b (cyan) at τ = 0 − leads to a complete loss of entanglement. Features like entanglement sudden death and revival observable in the subsequent dynamics of EN (for symmetric decay) can be traced back to the presence of the two-level system and vanish for strong asymmetry (γa/γ b = 0.3).
Turning toward asymmetric decay rates, γ a /γ b = 0.3, the degree of bipartite entanglement is generally lowered since the two-level system is much less populated. Furthermore, the delayed revival of entanglement is not that pronounced anymore as the probability that the two-level system is still excited after an emission event from the sixteen-level system is reduced. The fact that the vanishing of entanglement for a finite time span is lifted here completely indicates that this phenomenon can be traced back to the harsh state space restriction in oscillator b, which has here a less dramatic effect due to the large γ b .

E. Map of entanglement
The ratio γ a /γ b of the decay rates, as seen in the previous section, has a considerable impact on the characteristics of entanglement. Once the dimension of effective Hilbert space is determined, i.e., both the number of active oscillators and the number of accessible levels in each of these oscillators is fixed, there still remain two parameters, E J /E c2 J and γ a /γ b , which determine whether entanglement is present or not. This can be visualized in an entanglement map dividing the two-dimensional space spanned by these parameters into regions of separable or entangled states.
In the two-qubit case, used in Fig. 5(a) as a simple example, the PPT criterion already provides a necessary and sufficient condition for entanglement. Immediately, the corresponding steady-state entanglement map in Fig. 5(a) results, with entangled states (NPT) below and separable states (PPT) above E J /E c2 J = (γ 2 a + γ 2 b )/(γ a γ b ). The phenomenon of such an entanglement-separability threshold has also been observed in a number of similar steady-state systems, see, e.g., Refs. [63][64][65][66][67].
Determining the exact boundary line between entangled and separable states for an arbitrary bipartite system is usually more difficult due to the lack of criteria which are both necessary and sufficient for the detection of entanglement. Nonetheless, the basic structure of such an entanglement map in the bipartite case is always quite trivial as any state falls in either of only two classes. This, however, changes drastically when turning toward multipartite systems in the following section, where the entanglement structure naturally becomes much more complex [cf. Fig. 5(b)].

IV. MULTIPARTITE ENTANGLEMENT
The structure of entanglement between more than two parties is much richer than in the bipartite case and does not represent a trivial extension of these results [1][2][3][4][5][6]. Besides the basic question whether a given multipartite state is separable or entangled, we now have to specify the type of entanglement. A system consisting of N parties in a nonseparable state does not necessarily have to be genuinely N -partite entangled but can also be mseparable (1 < m < N ), i.e., there exists a splitting of the N parties into m groupings which are separable from each other. In addition to this, even genuinely N -partite entangled states can further be divided into different subclasses. The classification of multipartite states, in particular multipartite mixed states, is still a matter of intense research and, despite all efforts, it is far from being completely understood in general and is only explicitly known for very special cases.
Here, we therefore concentrate on one of these special cases, namely, a three-qubit system realized for κ q = 2 in our setup. In Sec. IV A, we briefly review a commonly used classification scheme for mixed three-qubit states before we apply this to our specific 2 × 2 × 2 system in Sec. IV B and chart the corresponding map of entanglement classes.

A. Classification of mixed three-qubit states
For pure three-qubit states, it is well known that there exist different equivalence classes of entanglement [9], and such a classification scheme can also be extended to mixed states [10].
A pure three-qubit state is called a fully separable state if it can be written as |ϕ fs = |ψ a ⊗ |ψ b ⊗ |ψ c . A biseparable state is a state where one of the three parties is separable from the other two. Consequently, there exist three different classes of biseparable states depending on which of the two subsystems are grouped together: for example, |ϕ bs,a = |ψ a ⊗ |ψ bc indicates a state where the parties b and c may be entangled. States which are  Figure 5. (a) Map of entanglement classes in the parameter space of the 2 × 2 system (κa = κ b = 2) in steady state spanned by the driving, EJ /E c2 J , and the ratio of decay rates, γa/γ b . Since the PPT criterion is both necessary and sufficient for the detection of entanglement here, any state can directly be assigned to either the class of entangled or separable states. (b) The corresponding map for the 2 × 2 × 2 system (κa = κ b = κc = 2) with one differing decay rate, γ a/b /γc = 1. Combining the statements of different entanglement criteria (C1-C4) reveals a rich structure of different regions in parameter space which are associated with one (or several) three-qubit entanglement classes.
neither fully separable nor biseparable are called genuinely tripartite entangled. These states, however, can further be divided into two classes of inequivalent states, the so-called W-class states |ϕ W with representative |W = (|001 + |010 + |100 )/ √ 3 and GHZ-class states |ϕ GHZ represented by |GHZ = (|000 + |111 )/ √ 2. Genuinely entangled states belonging to the GHZ class or W class cannot be transformed into one another by stochastic local operations and classical communication (SLOCC).
For mixed three-qubit states corresponding classes based on the results for pure states via convex combinations can be defined. A mixed state is considered as a fully separable state if it can be written as a convex combination of fully separable pure states |ϕ j fs , i.e., ρ fs = j p j |ϕ j fs ϕ j fs | with convex weights p j > 0 and j p j = 1. Accordingly, a mixed state is called biseparable if it can be expressed as a convex sum of pure biseparable states |ϕ j bs and pure fully separable states |ϕ j fs .
The pure biseparable states in this sum might be entangled with respect to different partitions. Additionally, we define subclasses of biseparable mixed states where pure biseparable states in the convex sum are all separable with respect to a fixed partition, i.e., either the partition a|bc, b|ac, or c|ab. States which are part of the W class can be expressed as a convex sum of pure W-class states |ϕ j W , pure biseparable states |ϕ j bs , and pure fully separable states |ϕ j fs . All mixed states which are not covered by the classification so far are in the GHZ class. States which belong either to the W or GHZ class are summarized as genuinely tripartite entangled states.

B. Map of entanglement classes
We will now discuss how to chart a map of entanglement classes for the various stationary states of our Josephson-photonics device in the three-qubit case, κ q = 2. Figure 5(b) shows the resulting map in the parameter space spanned by the driving strength and an additional parameter allowing for unequal damping of the cavities. To chart this map of entanglement classes, we will not directly rely on the above introduced definitions since a straightforward rewriting of the density matrix is usually not feasible. Instead, various witnesses or other criteria which have been developed based on these definitions are employed delivering statements on the entanglement nature of the stationary state for a certain region of parameter space. These multiple statements are then eventually combined to construct the full map.
We first find that the steady-state density matrix in the three-qubit case has a simple structure, when written in a standard product basis {|000 , |001 , |010 , . . . , |111 }. Alluding to the pattern of nonzero matrix elements, states of this structure are known as X states (here, there are some additional zeros). X states have been widely studied [68][69][70][71][72][73][74] with respect to entanglement and other quantum properties, in particular a subset of these states, which is called GHZ-diagonal [75]. Despite their simple structure, states of this form yield a rich pattern in the map of entanglement classes as we will see in the following.
Allowing for the decay rate of one of the qubits to differ from the other two, i.e., r = γ a/b /γ c = 1, the density matrix elements of ρ st in Eq. (6) can be explicitly found as simple functions of driving, With those expressions, we can now proceed in our entanglement classification focusing first on symmetrically damped oscillators, r = 1.
Here, we find the inequality violated for all states where E J /E c3 J < √ 2, implying genuine tripartite entanglement below this driving strength, see Fig. 5(b). • Criterion C2 based by Novo et al. [76] on a combination of PPT mixtures and permutationally invariant states, in contrast, detects states which are definitely not genuinely entangled but biseparable at the most. A PPT mixture is defined as a convex combination of PPT states with respect to a specific partition, i.e., ρ mix = p a ρ mix a +p b ρ mix b +p c ρ mix c with (ρ mix j ) Tj ≥ 0 ∀j. A permutationally invariant three-qubit state is biseparable if and only if it is a PPT mixture. Here, we can show that for √ 2 ≤ E J /E c3 J all states can be written as a convex sum of fully separable states and a permutationally invariant PPT mixture, i.e., all these states are not genuinely tripartite entangled. For √ 2 > E J /E c3 J , we can not make any statements.
• Criterion C3 follows work by Wölk et al. [45] to develop an entanglement witness based on the Hölder inequality which detects nonseparable states: If ρ is a fully separable state, then the inequality |ρ 000,111 | ≤ 4 √ ρ 000,000 ρ 011,011 ρ 101,101 ρ 110,110 holds. Violation is consequently a sufficient condition for ρ to be not fully separable, which indeed is the case here for all states where 62. This implies that ρ st is biseparable, but not fully separable, in the regime • Criterion C4 is one of a large number of powerful criteria that is available for states which are diagonal in the GHZ basis |Ψ ± j = (|ϕ a ϕ b ϕ c ± |φ aφbφc )/ √ 2 with ϕ l ∈ {0, 1} and ϕ l =φ l . Dür et al. [75] presented a necessary and sufficient criterion to decide whether such a GHZ-diagonal state is fully separable. We first apply an invertible local transformation to ρ st such that ρ st 000,111 = ρ st 111,000 ∈ R, thereby not changing the separability properties. For 6 ≤ E J /E c3 J , it can then be shown that all these states can be written as a convex sum of obviously fully separable states and GHZ-diagonal states which turn out to be fully separable. Again, we can make no statement for 6 > E J /E c3 J . For states in the regime 4.62 E J /E c3 J < 6, we can thus only state that they are not genuinely tripartite entangled.
Dropping the restriction to equally damped oscillators, for r = 1 the structure of entanglement becomes slightly modified as new regimes occur. Firstly, a gap arises between the regime boundaries based on C1 and C2 so that a region develops where states are known (by C3) to be not fully separable, but the type of entanglement could not be determined. Secondly, we can identify two additional regimes of biseparable states where the states are separable with respect to a fixed partition. This is shown again by rewriting ρ st as a sum of GHZ-diagonal states and fully separable states. Building on work of Dür et al. [75], we can prove that these GHZ-diagonal states are biseparable with respect to the corresponding partition. As apparent in Fig. 5(b), the overlap of these two regimes eventually forms the regime of fully separable states.
Intuitively one might think that a strong symmetry between the three oscillators would favor genuine tripartite entanglement. The results in Fig. 5(b) reveal, however, that in the near-symmetric case, r ≈ 1, the transition from genuine tripartite entanglement to biseparability and eventually further to full separability already occurs at comparatively moderate driving. Tuning the Josephson energy E J in this near-symmetric case allows to access a three-qubit state with the desired entanglement properties. For increasing asymmetry, the boundary lines between these regimes are shifted towards higher values of E J and we end up in an entangled state even for very strong driving. The blockade of further excitation and de-excitation processes once a photon has left the system from triple occupation |111 is particularly effective for strong asymmetry and protects the genuine tripartite entanglement against the impact of strong driving. The influence of the driving strength on the entanglement properties can again be understood by considering the relative size of the coherences |ρ st 000,111 | = |ρ st 111,000 | and the populations ρ st j,j (cf. explicit discussion for the 2 × 2 realization in Sec. III B). In the weak-driving limit, the coherences are dominant compared to the populations implying a strong quantum character in form of tripartite entanglement. Conversely, the coherences become small for strong driving, while the populations saturate leading to a more classical behavior which is reflected in full separability.
One important question which has not been addressed so far is whether the genuinely tripartite entangled states belong to the W or GHZ class. The simultaneous creation/annihilation processes of three photons by a single Cooper pair at the resonance considered obviously suggest that these states are part of the GHZ class. To show that a state lies within the GHZ class can in principle be proven by the use of witnesses; however, different commonly used witnesses [5,10,77] which we have tested on our specific 2 × 2 × 2 configuration have not detected GHZ-class states. Possibilities to show that the genuinely tripartite entangled states belong, contrary to our expectations, to the W class are unfortunately scarcely available. In particular, W-class states cannot be detected by witnesses in general as these are not designed to prove that a state lies inside a convex set [5].

V. EXPERIMENTAL SITUATION
Current experiments in Josephson photonics have already realized two-cavity setups and taken first steps in tackling the nonclassicality of microwave emission [26,27]. We now shortly turn to a discussion to what extent the promised versatility of Fig. 1(c) can actually be fulfilled and what improvements and modifications may still be necessary.
Constructing setups of three or more cavities and tuning the bias voltage V between the corresponding resonances to switch directly from a bipartite to a multipartite system does not seem to pose a principle challenge. Besides the number of active parties, also the Josephson energy E J , i.e., the class of entanglement, is easily tunable over a few orders of magnitude via the magnetic flux when using a SQUID configuration for the Josephson junction.
By contrast, the dimension of Hilbert space determined by the parameters κ q is currently fixed by design. While earlier experimental realizations were limited to the lowimpedance regime κ q 1, recent progress already al-  Figure 6. Witnesses Wj st = (a † ) j a j st − | (a † b † ) j st| for j = 1, 2, and 3 detecting steady-state entanglement in the symmetric 16 × 16 system (κa = κ b ≈ 0.23) for values of EJ /E c2 J where Wj st < 0. All three witnesses fail in detecting nonseparability in a regime of moderate driving between EJ /E c2 J ≈ 1.96 and 5.97, where the presence of entanglement has already been proven by the PPT criterion (cf. Fig. 2).
lows to reach values up to κ q ≈ 1.6 [78]. Achieving a κ q value which exactly matches one of the special values restricting the Hilbert space to a finite number of levels is actually not that crucial. Assuming N × M systems above is convenient for calculations and interpretation. However, except for extremely strong driving there will hardly be a difference between a near suppression and an exact vanishing of a certain transition matrix element [cf. Fig. 2(a) and also the discussion in Ref. [39], showing that deviations in the populations in a 2 × 2 system are of order O(δκ 2 )]. Nonetheless, if tunability is desired, employing SQUIDs at the end of a cavity or meta-material striplines constructed from long SQUID arrays [79], the effective cavity length (and therefore the resonance frequency) can be changed by magnetic fluxes. Moreover, the actual defining equation , where E C = 2e 2 /C q is the charging energy, corresponds to κ q = (2e 2 / ) L q /C q for the fundamental λ/4 mode ω (0) q . Accessing higher modes instead, therefore also gives immediate in-situ access to lower κ q values within currently existing devices.
Our theoretical investigations have been restricted to a small portion of the "space of entanglement phenomena" [Fig. 1(c)], namely, the full bipartite case and the low-dimensional 2 × 2 × 2 Hilbert space in the tripartite case, mainly due to the fact that only here a classification scheme is known and implementable criteria to detect these classes are available. The criteria used in our studies of bipartite and multipartite entanglement in the previous two sections, moreover, all require the full knowledge of the density matrix ρ. In an experiment, ρ can be reconstructed on the basis of quantum state tomography [80], where the density matrix of an unknown state is fully determined by repeatedly performing measurements  in different bases on an ensemble of identical copies of this state. In principle, investigating steady-state entanglement, here, such identical copies come for free. However, the whole procedure is challenging and, moreover, costly in terms of time, which as will be discussed below will pose severe limitations on implementing it in our setup. Particularly for systems with large Hilbert spaces, the knowledge of a quantum state ρ in an experimental situation will therefore usually be incomplete as based on a small number of observables only. For that reason, a large number of entanglement criteria in terms of directly measurable observables have been proposed. Sometimes these are simple witnesses relying on a single inequality only (see, e.g., Refs. [44,45,71,[81][82][83][84]), but also more elaborate detection schemes are known which, e.g., relate the PPT criterion of a given state ρ to the positivity of a corresponding (infinite) matrix of moments [85,86]. Whether one of these criteria is successful in detecting entanglement or not strongly depends on the structure of ρ. In principle, however, there exists an entanglement witness for any entangled state ρ [47]. To illustrate the limitations of simple witnesses in detecting entanglement, we apply in Fig. 6 the first three witnesses, j ∈ {1, 2, 3}, of the family [44] to the steady state ρ st of the symmetric 16 × 16 system. They detect entanglement, W j st < 0, in the weak-driving and strong-driving regime only, however, not in an intermediate regime between E J /E c2 J ≈ 1.96 and 5.97, where entanglement indeed is present (cf. Fig. 2).
Leakage of excited photons from the resonators into the electromagnetic environment is the dominant but not the only source of decoherence in an actual experimental realization. As discussed in Refs. [28], local voltage fluctuations at the Josephson junction associated with a rate γ J are comparatively weak, γ J /γ q 1, in experiments [21,24]. While they can be safely disregarded for some observables [28,32,39], they have a non-negligible impact on the entanglement properties in steady state. Describing the effect of voltage noise requires an extended model [28,32] with an extra degree of freedom for the number of Cooper pairs N that have passed the tunnel element. The bracketed term in the Hamiltonian [Eq. (1)] is here replaced by ( q e iη a † q + q e −iη a q ) with e iη = N |N N + 1|, where the phase difference η across the junction and N form a pair of conjugated variables, [η, N ] = i. The quantum master equation [Eq. (2)] is then modified by the additional dissipator γ J (2N ρN − N 2 ρ − ρN 2 )/2, including the decohering effects of voltage fluctuations on the system's dynamics. To analyze how voltage noise affects the entanglement properties in steady state, we study in Fig. 7 the logarithmic negativity E N for a 3 × 3 system and symmetric decay rates (γ 1 = γ 2 = γ) with γ J /γ = 0.01 during a switching-on/off process of the Josephson-photonics device. We assume for ρ(τ = 0) = ρ 0 that the two cavities are initially in their ground state and there is a welldefined phase difference across the Josephson junction, (e ±iη ) j = 1 with j ∈ N. Figure 7(a) pictures E N as a function of time for different values of the driving strength E J after switching on the device at τ = 0 by instantaneously turning up the driving from E J /E c2 J = 0 to the respective value. In practice, the switching process can be performed considerably quicker than other relevant timescales, in particular 1/γ [37].
The logarithmic negativity increases immediately and approaches, on a time scale of a few 1/γ, a quasistationary state decaying exponentially but only very slowly with rates of the order γ J . At γτ = 50, we finally switch of the device again by instantaneously reducing the driving to zero, resulting in a fast vanishing of the entanglement on time scales associated with 1/γ. In Fig. 7(b), E N is additionally plotted as a function of the driving E J /E c2 J at different moments in time after switching on the device. At γτ ≈ 4.0, the system has already reached its quasistationary state and the corresponding logarithmic negativity almost coincides with the reference values for vanishing γ J , where the system is in a true stationary state (dashed line, cf. Fig. 2). Due to the impact of voltage noise, the degree of entanglement is henceforth continuously reduced and after a long, but finite, waiting time the system would eventually end up in a completely disentangled state. Nonetheless, we can conclude that for weak voltage noise there exists a time span sufficiently long for experimental observation where the multifaceted entanglement properties discussed above for vanishing γ J are definitely present before switching off and restarting the device becomes necessary.

VI. CONCLUSIONS AND OUTLOOK
To study the complexity of entanglement phenomena in multipartite systems is a challenging and in many respects a still-unsolved problem, both in terms of theoretical concepts and in terms of experimental realizations.
In this paper, we demonstrated that Josephson-photonics devices which combine the basic elements of circuit QED in form of a voltage-biased Josephson junction interacting with an array of microwave cavities may serve as simple and versatile sources for entangled microwave photons. The way these photonic states are created is strikingly simple and does not require complicated pulse shaping. Due to their steady-state nature, they are relatively robust against perturbations detrimental to entanglement. Thus, this work sets the theoretical ground for Josephson photonics setups as testbeds for advanced classifications of bi-and multipartite entanglement. The versatility of these devices is offered by different experimental knobs which allow to change the structure of the effective Hamiltonian and the value of its parameters so that the system reaches widely diverse stationary states. We showed how to vary experimentally accessible parameters so that the number of entangled parties and the structure of Hilbert space for each parties can be chosen and entangled states belonging to different classes can be realized. In that manner, Josephson-photonics devices may be used to map out a considerable portion of the space of entangled states by changing one or several of these parameters as was demonstrated for the bi-and multipartite qubit case. The rich and complex entanglement phenomena encountered here manifest the necessity to employ advanced entanglement criteria developed in abstract entanglement theory in recent years. This gives support to the notion that the entanglement character of states naturally resulting from a simple driven, dissipa-tive dynamics may constitute a fruitful field for further research.
For future experiments, one important question will be the best choice of witness. Complexity or ease of measurement will have to be balanced with the power of various possible witnesses to distinguish between different classes and detect entanglement in a wider or more restricted range. This links to the challenge of mitigating the debilitating impact of low-frequency noise on entanglement; either by choosing witnesses which can be measured during the time of quasistationarity, by further improving phase stability possibly via phaselocking schemes, or by exploiting observables and measurement schemes insensitive to phase noise (cf. Fransoninterferometric schemes [87,88]).
In that context, we want to emphasize that in the current work we concentrated on questions concerning the entanglement of cavity modes. Various closely related questions can be studied in the frequency-and/or timedomain directly for output modes, including detailed modeling of filtering and photon detection. Related to these issues are possible limitations of our results inherent to the modeling of dissipation by a standard quantumoptical Lindbladian. How and when to improve on this modeling is an important topic for further studies.
For applications as an entanglement source, an important feature distinguishing Josephson-photonics devices from other circuit-QED setups is the continuous mode of operation, which suggests particular potential as highintensity source. Operating several Josephson junctions in a parallel, self-synchronized manner may conceivable serve this purpose.
For abstract entanglement theory, Josephson photonics constitute one powerful example that even simple Hamiltonians can dynamically generate a wealth of complex entanglement phenomena, worthy of interest. This may trigger efforts to characterize and understand the entanglement properties of states which are naturally occurring (as steady states) of similar dynamical systems as studied here; thereby complementing those where the investigated states are chosen arbitrarily or according to other criteria. Particularly fascinating is the question to what extent it is possible to directly connect the structure of the generating Hamiltonian to the entanglement class of the eventually resulting mixed states. For instance, instead of a tripartite Hamiltonian of the a † b † c † + c.c type as at the 2eV = (ω a + ω b + ω c ) resonance of Josephsonphotonics, we may consider a a † b † + b † c † + a † c † + c.c Hamiltonian, which is not realizable in pure form in the current setup. [There would be competing (a † ) 2 terms among other corrections.] We may then ask whether this can be related to entanglement of GHZ or W type, respectively, and try to pose and answer similar questions for other multipartite systems.