Superconducting Josephson-Based Metamaterials for Quantum-Limited Parametric Amplification: A Review

Abstract

In the last few years, several groups have proposed and developed their own platforms demonstrating quantum-limited linear parametric amplification, with evident applications in quantum information and computation, electrical and optical metrology, radio astronomy, and basic physics concerning axion detection. Here, we propose a short review on the physics behind parametric amplification via metamaterials composed by coplanar waveguides embedding several Josephson junctions. We present and compare different schemes that exploit the nonlinearity of the Josephson current-phase relation to mix the so-called signal, idler, and pump tones. The chapter then presents and compares three different theoretical models, developed in the last few years, to predict the dynamics of these nonlinear systems in the particular case of a 4-wave mixing process and under the degenerate undepleted pump assumption. We will demonstrate that, under the same assumption, all the results are comparable in terms of amplification of the output fields.

1. Introduction

In the last decade, microwave quantum electronics received a substantial boost by the advancements in superconducting circuits and dilution refrigerators technologies. These platforms allow experiments to be easily carried out in the mK regime, where the detection and manipulation of signals in the range $3-12$ GHz reaches energy sensitivities comparable to a single photon [1].

Solid state microwave quantum electronics is founded on a building block that has no analogous in quantum optics: the Josephson junction [2]. This, in fact, is a unique nondissipative and nonlinear component that represents the key element of a large series of quantum experiments.

Furthermore, microwave quantum electronics allows the exploration of the so-called ultrastrong coupling regime [3], hard to be reached in quantum optics, and it is worth mentioning that nonlinear resonator can be exploited to access relativistic quantum effects and quantum vacuum effects. To give an example, the Lamb shift [4] effect has been observed in superconducting artificial atom [5], while the dynamical Casimir effect [6, 7] has been promoted by properly engineered superconducting waveguide [8].

From the very beginning, superconducting electronics has been pushed by the strong interest coming from the quantum computation and information community. However, it has been only recently shown that a new concept of 1D metamaterial with embedded several Josephson junctions enables strong photon-photon on-chip interactions [9], allowing experimentalists to engineer dispersion relations that drive the waves traveling along artificial waveguides [10, 11]. These concepts and technologies allow the control and tunability of the wave mixing process. As an example, a weak signal traveling in a metamaterial can interact with a strong pump tone at a different frequency, activating the so-called parametric amplification [12]. The class of devices where these phenomena are promoted is commonly known as traveling-wave Josephson parametric amplifiers (TWJPA) and represents the solid state analogous to optical ${\chi }^n$ nonlinear crystals [13].

It has been shown that TWJPAs can act as quantum parametric amplifiers by reaching the so-called quantum limit [14]. With the purpose of a comparison to the state-of-the-art commercially available low-noise amplifiers, these latter can operate at $\omega /2\pi =4\mathrm{GHz}$ adding $k_B{T}_n/\mathrm{\hslash }\omega \approx 10$ noise photons having a noise temperature of $T_n=2\mathrm{K}$, while Josephson-based amplifiers can reduce this added noise up to $1/2$ photon, or even 0, depending on its working configuration.

The capability to beat the quantum limit is related to the so-called phase-sensitive amplification process, where the metamaterial can operate in degenerate mode (degenerate parametric amplifier, DPA), acting on two waves (signal and idler) at the same frequency (${\omega }_s={\omega }_i$) by amplifying and de-amplifying their position and momentum quadratures, respectively. In this view, DPA enables the preparation of squeezed states in the microwave regime. Even in the nondegenerate mode (nondegenerate parametric amplifier, NDPA, i.e., ${\omega }_s\ne {\omega }_i$), the phase-preserving nature of the quantum parametric amplification results in the entanglement condition among the signal and idler generated photons, composing a two-mode squeezed state [15]. It is worth mentioning how such a quantum state is an example of Einstein-Podolsky-Rosen state [16], where correlations between signal and idler are stronger than that allowed by classical theory [17].

It should be evident how superconducting electronics not only has demonstrated to be an ideal platform for microwave quantum parametric amplification but also has pushed forward the research field focusing on the generation of nonclassical radiation with attractive potential applications in metrology and quantum information processing.

2. Historical evolution of the traveling-wave parametric amplifiers

The theory of a new concept of microwave amplifier was developed by Cullen [12] in 1959. In his paper, Cullen showed a novel mechanism of periodic transfer of power between a pump tone and a signal traveling in a transmission line composed of a voltage dependent capacitance per unit length. A nonlinear component of an RLC circuit can change periodically the resonance frequency of the whole system, leading to a novel way of making broadband amplification, the so-called parametric amplification. In Figure 1, we report two toy models for parametric amplification in mechanical systems with their electrical counterparts.

Figure 1.

(a) Sketch of a swing process. An oscillating system at a frequency ${\omega }_s$ is excited by parametric amplification via periodical changes of the center of mass position at a frequency ${\omega }_p=2{\omega }_s$. (b) LC circuit with variable (nonlinear) C and L components. The case in which the capacitance C is periodically changed in time is the circuit analogous to the mechanical system represented in (a), while the case having an oscillating inductance L mimics the condition sketched in (c), consisting in a torque pendulum with variable inertia momentum [18].

One of the first realizations of Cullen’s idea was made by Mavaddat et al. in 1962 [19]. The signal line was basically a low pass filter, in which the shunt elements were similar varactor diodes. There, the nonlinearity was given by the specific capacitance-voltage relation of the varactor diodes, which is highly nonlinear for relative small voltage values. In this pioneering experiment, a gain of 10 dB and a bandwidth of 3 MHz were shown.

After the theorizing and the subsequent discovery of the Josephson effect [2], it was understood that an easy way to embed a nonlinear component into a transmission line and simultaneously reduce losses was to build a nonlinear inductance made of superconducting material, exploiting a Josephson junction as a source of nonlinearity following the vanguard idea by Sweeny and Mahler [13]. There, the parametric amplifier was modeled by a first-order small-signal theory with the same approach adopted to predict the behavior of GaAsFET transmission line amplifiers. The proposed design consisted of a superconducting thin-film niobium transmission line, composed by a coplanar waveguide integrating a large number of Josephson junctions.

The first realization of a traveling-wave parametric amplifier embedding a series of Josephson junctions was possible due to the PARTS process developed at IBM [9]. Exploiting niobium/aluminum technology, Yurke et al. [20] reported the construction and characterization of a coplanar waveguide, in which the central trace was composed by an array of 1000 Josephson junctions. The experiment was there performed in reflection mode, by terminating one end of the device with a short, leading to a relative high gain of 16 dB but a narrow bandwidth of 125 MHz and a noise temperature of 0.5 ± 0.1 K. The mismatch between the theoretical model and the experimental data has resulted in the understanding of a lack of a complete description of the physics behind this device when working in a small-signal regime. The study of the collective behavior of groups of Josephson junctions forming a transmission line has been an active field of study of several theoretical works [21, 22]. Subsequently, the use of numerical analysis [23] helped in clarifying how wave propagation acts inside this kind of transmission line, giving information on cutoff propagation, dispersive behavior, and shock-wave formation. An analytical model of a Josephson traveling-wave amplifier of greater complexity was developed by Yaakobi et al. [24]. There, a transmission line made of a series of capacitively shunted Josephson junctions was considered.

One of the main limitations concerning the maximum achievable gain, common to all the TWJPAs concepts, is represented by the phase mismatch between the different tones into the line. In particular, it has to be noticed that even though the incoming waves can be in phase, photon-photon interactions between different tones (cross-phase modulation) or the same tone (self-phase modulation) lead to a modification of the phase of the traveling tones themselves. Indeed, quantum mechanically speaking, the power transport between the pump and the signal waves takes place through a photon energy conversion between the pump and the signal. This means that for an efficient energy exchange, conservation of both energy and momentum needs to take place. The latter condition is the corpuscular analogous to the phase matching requirement between the different electromagnetic waves. An engineering solution to overcome this problem is represented by the so-called resonant phase matching (RPM) [11]. O’Brien et al. analyzed this method theoretically on a simple transmission line made of a series of Josephson junctions capacitively shunted to ground operating in the so-called 4-wave mixing (4WM) regime. In their model, they shunted the transmission line with several LC resonators with a resonance frequency slightly above the pump tone. Doing this, they were able to show the rise of a stop band in the dispersion relation, which is able to re-phase the pump with the signal tones by changing the pump wave vector, favoring the wave mixing.

O’Brien’s design was realized not long after [10] using Al technology. In their design, the unit cell of the transmission line was composed by three single nonlinear Josephson cells, the shunt capacitor was made using low-loss amorphous silicon dielectric and a resonator was placed after each group of 17 unit cells. The device showed a maximum gain of 12 dB over a 4 GHz bandwidth centered on $\approx 5\mathrm{GHz}$. Moreover, the authors explain that variations of 2–3 dB in the gain most likely come from imperfect impedance matching between sections and at the level of the bond pads.

A similar design was adopted by Macklin et al. [14] to prove experimentally the capability of a TWJPA combined with the RPM technique to be used as a reliable tool for qubits readout. In this paper, the TWJPA, based on Nb technology and a different RPM periodicity, was first characterized, showing a gain of 20 dB over a 3 GHz bandwidth. Moreover, the quantum efficiency of the amplifier was tested when coupled with a 3D transmon qubit, leading to an efficiency value of $0.49\pm 0.01$. A key point of this experiment was the proof that a single TWJPA could be able to perform the readout of more than 20 qubits, thanks to its high dynamic range and multiplexing capabilities. RPM has shown remarkable capabilities and is a promising technique to overcome phase mismatch. It can be implemented in multiple ways [25], by the way, it has to be noticed that this method requires an increase of design complexity, lower tolerances on the constructing parameters, and longer propagation lengths (2 cm–1 m).

Another option to solve the mismatch problem was suggested by Bell and Samlov [26], who proposed a self-phase matching transmission line embedding a series of asymmetric superconducting quantum interference devices (SQUIDs). The remarkable feature of this design is that it does not need any resonant circuit to achieve phase matching. This TWJPA is indeed able to tune the nonlinearity of its SQUIDs just through the use of an external magnetic field. Zhang et al. realized that this design [27] proves the wide tunability on positive and negative values of the Kerr nonlinearity by a magnetic flux and its capability to assist phase matching in the 4WM process. The 4WM process is intrinsically affected by phase mismatch because it takes origin from a cubic (Kerr-like) nonlinearity of the current-phase relation of the SQUIDs composing the TWJPA, getting unwanted effects from self-phase and cross-phase modulations.

Zorin showed [28] that by embedding a chain of rf-SQUIDs into a coplanar waveguide, it is possible to tune both the second and third order nonlinearities of their phase-current relation. This is a totally a novel approach to the TWJPA, since the possibility to use a quadratic term as a source of nonlinearity, allows to work in the 3-Wave Mixing regime (3WM), as theorized by Cullen 57 years before. It is well-known that 3WM has several advantages when compared to 4WM. Firstly, it allows to operate with a minimal phase mismatch. Secondly, it requires a lesser pump power to achieve the same amplification per unit length. Eventually, it separates signal and idler from pump tones, easing the engineering of the experimental setup by removing the requirement of heavy filtering in the middle of the amplification band. A proof of principle based on the Zorin’s layout [29] showed a gain reaching 11 dB over a 3 GHz bandwidth.

A step forward in controlling the metamaterial nonlinearities was attempted by Miano et al. [30] achieving an independent tune of both second and third order terms in the current-phase relation by adjusting the bias current in some inductive circuits surrounding the transmission line. This technology takes the name of symmetric traveling-wave parametric amplifier (STWPA), its peculiarity arising from the symmetric arrangement of the rf-SQUIDs that compose the transmission line. This device concept represents the state-of-the-art in the field, allowing the exploration a wide portion of the control parameters space, leading to a maximum estimated gain of 17 dB and a 4 GHz bandwidth.

3. Theoretical models for a 4WM process in a TWJPA

In the last decade, different theoretical models have been developed to predict the behavior of an electric transmission line containing an array of Josephson junctions, employed as nonlinear elements. In this section, we will focus on those models developed to predict the behavior of a TWJPA in the particular case of a 4WM process, under undepleted degenerate pump approximation (i.e., assuming that the power held by the pump wave is at first approximation constant and larger than the one owned by the signal and the idler). We will firstly focus on the classical theory proposed by Yaakobi et al. in 2013 [24] and O’Brien et al. in 2014 [11], in which the behavior of the transmission line is derived imposing the current conservation in the system. This starting assumption leads to the definition of a partial differential nonlinear equation that can be turned into a system of coupled mode equations, providing the expression of the amplitude of the pump, signal, and idler tones along the transmission line. Subsequently, we will discuss two different quantum approaches for the description of the parametric amplifier dynamics. The first one, proposed by Grimsmo and Blais in 2017 [31], exploits a Hamiltonian based on continuous-mode operators to derive, in an interaction picture frame, a device’s output field. The second one instead, proposed by van der Reep in 2019 [32], derives a system of coupled mode equations for the creation and annihilation quantum operators starting from a Hamiltonian based on discrete-mode operators.

The theories presented in this chapter will be based on a series of simplifying assumptions, whose experimental realization could be difficult to be obtained. For instance, in a real device, the undepleted pump approximation is hardly respected along the entire extension of the device because, along the line, the pump tone transfers a non-negligible amount of energy to the signal and idler one. The depletion effects, resulting in a reduction of the gain and of the dynamics-range of the amplifier, have been studied both in a classical and quantum frame [28, 33].

In all of these models, a lossless electrical circuit composed by the repetition of an elementary cell, whose structure is shown in Figure 2, is taken into account. In order to standardize the notations, we assume that the Josephson junctions embedded in the transmission line are identical (i.e., they have the same critical current $I_c$) and that the current flowing through the n-th junction can be expressed through the nonlinear relation

Figure 2.

Electrical equivalent representation of a repetition of Josephson junctions embedded in a transmission line. The junctions are modeled as an LC resonant circuit. The length of the unitary cell length is represented by a.

$$I_{J,n}=I_{c}sin\left({\phi }_n\right)$$E1

where ${\phi }_n\left(t\right)$ is the phase difference across the junction of the macroscopic wave functions of the two superconductive electrodes. The relation between ${\phi }_n$ and the voltage-difference across the junction is given by the Faraday’s induction law

$$\mathrm{\Delta }{V}_n=V_{n+1}-V_n=\frac{{\mathrm{\Phi }}_0}{2\pi }\frac{d{\phi }_n}{\mathit{dt}}=\frac{{\mathrm{\Phi }}_0}{2\pi }\frac{d}{\mathit{dt}}\left[{\varphi }_{n+1}-{\varphi }_n\right]$$E2

where ${\mathrm{\Phi }}_0=h/\left(2e\right)$ is the magnetic quantum flux (with h the Planck constant and e the elementary charge), whereas ${\varphi }_n\left(t\right)$ is the absolute phase in the n-th node of the circuit. The phase at the n-th node (${\varphi }_n$) can be converted into a flux at the n-th node, and vice versa, through the relation ${\mathrm{\Phi }}_n=\left({\mathrm{\Phi }}_0/2\pi \right){\varphi }_n$.

Furthermore, we define $C_J$ the capacitance associated to the n-th Josephson junction and $L_{J,n}$ its inductance, defined as

$$L_{J,n}=\frac{\mathrm{\Delta }{V}_n}{{\mathit{dI}}_{J,n}/\mathit{dt}}=\frac{{\mathrm{\Phi }}_0}{2\pi }\frac{1}{I_{c}cos\left({\phi }_n\right)}=\frac{L_{J_0}}{cos\left({\phi }_n\right)}$$E3

where $L_{J_0}={\mathrm{\Phi }}_0/\left(2\pi I_c\right)$ is the inductance of the Josephson junction for a phase difference ${\phi }_n=0$.

The energy stored in the n-th Josephson junction can be expressed, using the definitions given in Eqs. (1) and (2), as

$$U_{J,n}={\int }_{t_0}^t\mathit{VI}d{t}^{\prime }={\int }_{t_0}^t\frac{{\mathrm{\Phi }}_0}{2\pi }\frac{d{\phi }_n\left(t^{\prime }\right)}{d{t}^{\prime }}{I}_{c}sin\left({\phi }_n\right)d{t}^{\prime }=I_c\frac{{\mathrm{\Phi }}_0}{2\pi }\left[1-\mathit{cos}\left({\phi }_n\left(t\right)\right)\right],$$E4

under the assumption that ${\phi }_n\left(t_0\right)=0$, and approximated through a first-order power expansion as

$$U_{J,n}=I_c\frac{{\mathrm{\Phi }}_0}{2\pi }\left[1-\mathit{cos}\left({\phi }_n\right)\right]=\frac{1}{2L_{J_0}}{\mathrm{\Delta \Phi }}_n^2-\frac{1}{24L_{J_0}}{\left(\frac{2\pi }{{\mathrm{\Phi }}_0}\right)}^2{\mathrm{\Delta \Phi }}_n^4+O\left({\mathrm{\Delta \Phi }}_n^6\right)$$E5

where ${\mathrm{\Delta \Phi }}_n={\mathrm{\Phi }}_{n+1}-{\mathrm{\Phi }}_n$.

Finally, we assume identical coupling capacitances $C_0$ between the transmission line and ground.

3.1 The classical theoretical model

In this subsection, we will present the main steps for the derivation of the classical model presented in [11, 24]. Under proper assumption, this model allows to determine analytically the amplitude of the signal’s and idler’s waves along the transmission line.

Expressing the current flowing through each branch of the circuit presented in Figure 2 in terms of absolute phases ${\varphi }_n$ and imposing the current conservation in the n-th node (i.e., $I_{J,n-1}+I_{C_J,n-1}=I_{J,n}+I_{C_J,n}+I_{C_0,n}$), a differential equation for the absolute phase, in the case of a weak nonlinearity, can be obtained:

$$\begin{array}{c}-C_0\frac{d^2}{{\mathit{dt}}^2}\left[{\varphi }_n\right]=-C_J\frac{d^2}{{\mathit{dt}}^2}\left[{\varphi }_{n+1}+{\varphi }_{n-1}-2{\varphi }_n\right]-\frac{1}{L_{J_0}}\left[{\varphi }_{n+1}+{\varphi }_{n-1}-2{\varphi }_n\right]\\ +\frac{{\mathrm{\Phi }}_0}{2\pi }\frac{1}{6I_c^2{L}_{J_0}^3}\left[{\left({\varphi }_{n+1}-{\varphi }_n\right)}^3-{\left({\varphi }_n-{\varphi }_{n-1}\right)}^3\right]\end{array}$$E6

where the last term derives from the first-order approximation of the nonlinear behavior of the Josephson’s inductance.

Assuming the length a of the elementary cell is much smaller than the wave lengths of the propagating waves $\lambda$ (i.e., $a/\lambda \ll 1$), the discrete index n can be replaced by a continuous position x along the line (i.e., ${\varphi }_n\left(t\right)\to \varphi \left(x,t\right)$) and the phase differences can be expressed, at the second order approximation, as:

$${\varphi }_{n+1}-{\varphi }_n\approx a\frac{\partial \varphi }{\partial x}+\frac{1}{2}{a}^2\frac{{\partial }^2\varphi }{\partial x^2}$$E7

$${\varphi }_n-{\varphi }_{n-1}\approx a\frac{\partial \varphi }{\partial x}-\frac{1}{2}{a}^2\frac{{\partial }^2\varphi }{\partial x^2}$$E8

In this way, it is possible to define a nonlinear differential equation for the continuous absolute phase $\varphi \left(x,t\right)$:

$$C_0\frac{{\partial }^2\varphi }{\partial t^2}-\frac{a^2}{L_{J_0}}\frac{{\partial }^2\varphi }{\partial x^2}-C_J{a}^2\frac{{\partial }^4\varphi }{\partial x^2\partial t^2}=-\frac{a^4}{2I_c^2{L}_{J_0}^3}\frac{{\partial }^2\varphi }{\partial x^2}{\left(\frac{\partial \varphi }{\partial x}\right)}^2$$E9

In the case of a weakly nonlinear medium, the dispersion law can be derived from Eq. (9), considering the left-hand side being equal to zero and imposing a plane-wave solution $\varphi \left(x,t\right)\propto e^{i\left(\mathit{kx}-\mathit{\omega t}\right)}$:

$$k\left(\omega \right)=\frac{\omega \sqrt{L_{J_0}{C}_0}}{a\sqrt{1-L_{J_0}{C}_J{\omega }^2}}$$E10

The solutions of Eq. (9) can be expressed, as shown by O’Brien et al. [11], in the form of a superposition of three waves (pump, signal, and idler) whose amplitudes are complex functions of the position along the line:

$$\varphi \left(x,t\right)=\sum _{n=p,s,i}\mathit{Re}\left[A_n\left(x\right)e^{i\left(k_{n}x-{\omega }_{n}t\right)}\right]=\frac{1}{2}\sum _{n=s,i,p}\left[A_n\left(x\right)e^{i\left(k_{n}x-{\omega }_{n}t\right)}+c.c\right]$$E11

The case of a 4WM process with a degenerate pump can be taken into account by imposing the frequency matching condition $2{\omega }_p={\omega }_s+{\omega }_i$. Replacing this particular solution in Eq. (9) and assuming that, along the line, the amplitudes are slowly varying (i.e., $\mid {\partial }^2{A}_n/\partial x^2\mid \ll k_n\mid \partial A_n/\partial x\mid \ll k_n^2\mid A_n\mid$) and that ${\left|A_s\right|}^2$ and ${\left|A_i\right|}^2$ are negligible (i.e., $\mid A_{s,i}\mid \ll \mid A_p\mid$, strong pump approximation), we obtain a system of three coupled differential equations for the amplitudes $A_n\left(x\right)$ that describe the energy exchange between the three waves along the line:

$$\frac{\partial A_p}{\partial x}=i{\vartheta }_p{\left|A_p\right|}^2{A}_p+2{\mathit{iX}}_p{A}_p^{\ast }{A}_s{A}_i{e}^{i\mathrm{\Delta }kx}$$E12

$$\frac{\partial A_{s\left(i\right)}}{\partial x}=i{\vartheta }_{s\left(i\right)}{\left|A_p\right|}^2{A}_{s\left(i\right)}+{\mathit{iX}}_{s\left(i\right)}{A}_p^2{A}_{i\left(s\right)}{e}^{i\mathrm{\Delta }kx}$$E13

where $\mathrm{\Delta }k=2k_p-k_s-k_i$ is the chromatic dispersion. The term ${\vartheta }_p$ is responsible for the self-phase modulation of the pump tone, while ${\vartheta }_{s\left(i\right)}$ is responsible for the cross-phase modulation between the pump tone and the signal or idler, respectively. These terms can be expressed as

$${\vartheta }_p=\frac{a^4{k}_p^5}{16C_0{I}_c^2{L}_{J_0}^3{\omega }_p^2}\text{and}{\vartheta }_{s\left(i\right)}=\frac{a^4{k}_p^2{k}_{s\left(i\right)}^3}{8C_0{I}_c^2{L}_{J_0}^3{\omega }_{s\left(i\right)}^2}$$E14

while the coupling constants $X_n$, depending on the circuit parameters, are defined as

$$X_p=\frac{a^4{k}_p^2{k}_s{k}_i\left(k_p-\mathrm{\Delta }k\right)}{16C_0{I}_c^2{L}_{J_0}^3{\omega }_p^2}\text{and}{X}_{s\left(i\right)}=\frac{a^4{k}_p^2{k}_s{k}_i\left(k_{s\left(i\right)}+\mathrm{\Delta }k\right)}{16C_0{I}_c^2{L}_{J_0}^3{\omega }_{s\left(i\right)}^2}$$E15

Expressing the complex amplitudes $A_n\left(x\right)$ in a the co-rotating frame

$$A_n\left(x\right)=A_{n_0}{e}^{i{\vartheta }_n{\left|A_{p_0}\right|}^{2}x}$$E16

it can be demonstrated that, working under the undepleted pump approximation $\mid A_p\left(x\right)\mid =A_{p_0}\gg \mid A_{s\left(i\right)}\left(x\right)\mid$, the amplitude of the signal can be expressed as

$$A_s\left(x\right)=\left[A_{s_0}\left(cosh\left(g_{1}x\right)-\frac{i{\mathrm{\Psi }}_1}{2g_1}sinh\left(g_{1}x\right)\right)+i\frac{X_{s\left(i\right)}{\left|A_{p_0}\right|}^2}{g_1}{A}_{s{\left(i\right)}_0}^{\ast }sinh\left(g_{1}x\right)\right]e^{i\frac{{\mathrm{\Psi }}_1}{2}x}$$E17

where ${\mathrm{\Psi }}_1=\mathrm{\Delta }k+\left(2{\vartheta }_p-{\vartheta }_s-{\vartheta }_i\right){\left|A_{p_0}\right|}^2=\mathrm{\Delta }k+\vartheta {\left|A_{p_0}\right|}^2$ is the total phase mismatch and $g_1$ is the exponential complex gain factor, defined as

$$g_1=\sqrt{X_s{X}_i^{\ast }{\left|A_{p_0}\right|}^4-{\left(\frac{{\mathrm{\Psi }}_1}{2}\right)}^2}$$E18

The total gain of an amplifier, composed by the repetition of $N$ elementary cells, can then be expressed as $G_s\left(\mathit{\text{aN}}\right)={\left|A_s\left(\mathit{\text{aN}}\right)/A_{s_0}\right|}^2$.

3.2 Quantum Hamiltonian model based on continuous-mode operators

A standard method to treat quantum superconducting circuits is represented by the lumped element approach [34]. In this latter, the Hamiltonian of the quantum circuit is straightforwardly derived from its classical counterpart by promoting fields to operators and properly imposing commutating relations. In this view, one can proceed by deriving the Lagrangian of a TWJPA composed by the repetition of N unitary cells, under first nonlinear order approximation, as

$$\mathcal{L}=\sum _{n=0}^{N-1}\left[\frac{C_0}{2}{\left(\frac{\partial {\mathrm{\Phi }}_n}{\partial t}\right)}^2+\frac{C_J}{2}{\left(\frac{\partial {\mathrm{\Delta \Phi }}_n}{\partial t}\right)}^2-E_{J_0}\left(1-cos\left(\frac{2\pi }{{\mathrm{\Phi }}_0}{\mathrm{\Delta \Phi }}_n\right)\right)\right]\approx$$E19

$$\approx \sum _{n=0}^{N-1}\left[\frac{C_0}{2}{\left(\frac{\partial {\mathrm{\Phi }}_n}{\partial t}\right)}^2+\frac{C_J}{2}{\left(\frac{\partial {\mathrm{\Delta \Phi }}_n}{\partial t}\right)}^2-\frac{1}{2L_{J_0}}{\mathrm{\Delta \Phi }}_n^2-\frac{1}{24L_{J_0}}{\left(\frac{2\pi }{{\mathrm{\Phi }}_0}\right)}^2{\mathrm{\Delta \Phi }}_n^4\right]$$E20

where $E_{J_0}=I_c{\mathrm{\Phi }}_0/2\pi =I_c{L}_{J_0}$. Under the assumption that $a/\lambda \ll 1$ it is possible, as performed in the previous subsection, to replace the discrete index $n$ with a continuous position $x$ along the line (i.e., ${\mathrm{\Phi }}_n\left(t\right)\to \mathrm{\Phi }\left(x,t\right)$) and approximate, at the first order, ${\mathrm{\Delta \Phi }}_n\to a\mathrm{\partial \Phi }\left(x,t\right)/\partial x$. Furthermore, extending the system via two lossless semi-infinite transmission lines (characterized by a constant distributed capacitance $c_0$ and a constant distributed inductance $l_0$), the Lagrangian can be expressed through a space integral extending from $x=-\infty$ to $x=+\infty$ as

$$\mathcal{L}\left[\mathrm{\Phi },\frac{\mathrm{\partial \Phi }}{\partial t}\right]=\frac{1}{2}{\int }_{-\infty }^{\infty }\left[c\left(x\right){\left(\frac{\mathrm{\partial \Phi }}{\partial t}\right)}^2+\frac{1}{{\omega }_J^2\left(x\right)l\left(x\right)}\left(\frac{{\partial }^2\mathrm{\Phi }}{\partial x\partial t}\right)-\frac{1}{l\left(x\right)}{\left(\frac{\mathrm{\partial \Phi }}{\partial x}\right)}^2+\gamma \left(x\right){\left(\frac{\mathrm{\partial \Phi }}{\partial x}\right)}^4\right]\mathit{dx}$$E21

where $c\left(x\right)$ and $l\left(x\right)$ are the distributed capacitance and inductance of the system, defined as

$$c\left(x\right)=\left\{\begin{array}{ll}{c}_0& x<0\\ C_J/a& 0<x<z\\ c_0& x>z\end{array}\right.\text{and}l\left(x\right)=\left\{\begin{array}{ll}{l}_0& x<0\\ L_{J_0}/a& 0<x<z\\ l_0& x>z\end{array}\right.$$E22

$${\omega }_J\left(x\right)=\left\{\begin{array}{ll}\infty & x<0\\ 1/\sqrt{L_{J_0}{C}_J}& 0<x<z\\ \infty & x>z\end{array}\right.\text{and}\gamma \left(x\right)=\left\{\begin{array}{ll}0& x<0\\ \left(a^3{E}_{J_0}/12\right){\left(2\pi /{\mathrm{\Phi }}_0\right)}^4& 0<x<z\\ 0& x>z\end{array}\right.$$E23

where $z$ is the length of the TWJPA, ${\omega }_J\left(x\right)$ is the junction’s plasma frequency, and $\gamma \left(x\right)$ is the term deriving from the nonlinearity of the junctions.

From Eq. (21), one can easily derive the Euler-Lagrange equation whose form, for $0<x<z$, is equal to Eq. (9), giving the same dispersion relation (Eq. (10)) under analogous assumptions. Instead, outside the nonlinear region, the wave vector turns out to be $k_{\omega }\left(x\right)=\sqrt{c_0{l}_0}{\omega }^2$.

The Hamiltonian of the system can be derived from the Lagrangian by taking into account $\pi \left(x,t\right)=\delta \mathcal{L}/\delta \left[\mathrm{\partial \Phi }/\partial t\right]$, the canonical momentum of the flux $\mathrm{\Phi }\left(x,t\right)$:

$$\begin{array}{l}H\left[\mathrm{\Phi },\pi \right]={\int }_{-\infty }^{\infty }\left[\pi \frac{\mathrm{\partial \Phi }}{\partial t}\right]\mathit{dx}-\mathcal{L}\\ =\frac{1}{2}{\int }_{-\infty }^{\infty }\left[c\left(x\right){\left(\frac{\mathrm{\partial \Phi }}{\mathrm{\partial \Phi }t}\right)}^2+\frac{1}{l\left(x\right)}{\left(\frac{\mathrm{\partial \Phi }}{\partial x}\right)}^2+\frac{1}{{\omega }_p^2\left(x\right)l\left(x\right)}{\left(\frac{{\partial }^2\mathrm{\Phi }}{\partial x\partial t}\right)}^2\right]\mathit{dx}-\frac{\gamma }{2}{\int }_0^z\left[{\left(\frac{\mathrm{\partial \Phi }}{\partial x}\right)}^4\right]\mathit{dx}\\ =H_0+H_1\end{array}$$E24

where the term $H_0$ represents the linear contributions to the energy of the system, while $H_1$ is the first-order nonlinear contribution.

This Hamiltonian can be converted to its quantum form promoting the field $\mathrm{\Phi }\left(x,t\right)$ to the quantum operator $\widehat{\mathrm{\Phi }}\left(x,t\right)$.

In direct analogy with Eq. (11), one can express the flux operator in terms of continuous-mode functions [34], such as ${\widehat{H}}_0$ is diagonal in the plane-waves unperturbed modes decomposition:

$$\widehat{\mathrm{\Phi }}\left(x,t\right)=\sum _{\nu =L,R}{\int }_0^{\infty }\left[\sqrt{\frac{\mathrm{\hslash }l\left(x\right)}{4\pi k_{\omega }\left(x\right)}}{\widehat{a}}_{\mathit{\nu \omega }}{e}^{i\left(\pm k_{\omega }\left(x\right)x-\mathit{\omega t}\right)}+\mathrm{H}.\mathrm{c}.\right]\mathit{d\omega }$$E25

where the subscript $R$ denotes a progressive wave, while $L$ denotes a regressive wave (i.e., ${\widehat{a}}_{\mathit{R\omega }}$ represents the annihilation operator of a right-moving field of frequency $\omega$). In [35], it is demonstrated that by replacing the definition given in Eq. (25) into the linear Hamiltonian, ${\widehat{H}}_0$ takes the form

$${\widehat{H}}_0=\sum _{\nu =R,L}{\int }_0^{\infty }\left[\mathrm{\hslash }\omega {\widehat{a}}_{\mathit{\nu \omega }}^{†}{\widehat{a}}_{\mathit{\nu \omega }}\right]\mathit{d\omega }$$E26

(where the zero-point energy, which does not influence the dynamics of the amplifier, has been omitted).

Using the expansion of $\widehat{\mathrm{\Phi }}\left(x,t\right)$ introduced above, under the hypothesis of a strong right-moving classical pump centered in ${\omega }_p$, and that the fields ${\widehat{a}}_{\mathit{\nu \omega }}$ are small except for frequencies closed to the pump frequency (i.e., replacing ${\widehat{a}}_{\mathit{\nu \omega }}$ with ${\widehat{a}}_{\mathit{\nu \omega }}+b\left(\omega \right)$, where $b\left(\omega \right)$ is a complex valued function centered in ${\omega }_p$), the nonlinear Hamiltonian ${\widehat{H}}_1$ can be expressed under strong pump approximation, at the first order in $b\left(\omega \right)$, as the sum of three different contributions:

$${\widehat{H}}_1={\widehat{H}}_{\mathit{CPM}}+{\widehat{H}}_{\mathit{SQ}}+H_{\mathit{SPM}}$$E27

In the expressions of these contributions, the fast rotating terms and the highly phase mismatched left-moving field have been neglected:

$${\widehat{H}}_{\mathit{CPM}}=-\frac{\mathrm{\hslash }}{2\pi }{\int }_0^{\infty }d{\omega }_{s}d{\omega }_{i}d{\mathrm{\Omega }}_{p}d{\mathrm{\Omega }}_{p^{\prime }}\sqrt{k_{{\omega }_s}{k}_{{\omega }_i}}{\beta }^{\ast }\left({\mathrm{\Omega }}_p\right)\beta \left({\mathrm{\Omega }}_{p^{\prime }}\right)\mathrm{\Upsilon }\left({\omega }_s,{\omega }_i,{\mathrm{\Omega }}_p,{\mathrm{\Omega }}_{p^{\prime }}\right){\widehat{a}}_{R{\omega }_s}^{†}{\widehat{a}}_{R{\omega }_i}+\mathrm{H}.\mathrm{c}$$E28

describes the cross-phase modulation,

$${\widehat{H}}_{\mathit{SQ}}=-\frac{\mathrm{\hslash }}{4\pi }{\int }_0^{\infty }d{\omega }_{s}d{\omega }_{i}d{\mathrm{\Omega }}_{p}d{\mathrm{\Omega }}_{p^{\prime }}\sqrt{k_{{\omega }_s}{k}_{{\omega }_i}}\beta \left({\mathrm{\Omega }}_p\right)\beta \left({\mathrm{\Omega }}_{p^{\prime }}\right)\mathrm{\Upsilon }\left({\omega }_s,{\mathrm{\Omega }}_p,{\omega }_i,{\mathrm{\Omega }}_{p^{\prime }}\right){\widehat{a}}_{R{\omega }_s}^{†}{\widehat{a}}_{R{\omega }_i}^{†}+\mathrm{H}.\mathrm{c}.$$E29

describes the broadband squeezing, and

$$H_{\mathit{SPM}}=-\frac{\mathrm{\hslash }}{4\pi }{\int }_0^{\infty }d{\omega }_{s}d{\omega }_{i}d{\mathrm{\Omega }}_{p}d{\mathrm{\Omega }}_{p^{\prime }}\sqrt{k_{{\omega }_s}k{\omega }_i}{\beta }^{\ast }\left({\mathrm{\Omega }}_p\right)\beta \left({\mathrm{\Omega }}_{p^{\prime }}\right)\mathrm{\Upsilon }\left({\omega }_s,{\omega }_i,{\mathrm{\Omega }}_p,{\mathrm{\Omega }}_{p^{\prime }}\right)b^{\ast }\left({\omega }_s\right)b\left({\omega }_i\right)+\mathrm{H}.\mathrm{c}.$$E30

describes the self-phase modulation. $\beta \left(\mathrm{\Omega }\right)$ is the dimensionless pump amplitude, proportional to the ratio between the pump current $I_J\left({\mathrm{\Omega }}_p\right)$ and the critical current of the junctions $I_c$:

$$\beta \left({\mathrm{\Omega }}_p\right)=\frac{I_J\left({\mathrm{\Omega }}_p\right)}{4I_c}$$E31

The function $\mathrm{\Upsilon }\left({\omega }_1,{\omega }_2,{\omega }_3,{\omega }_4\right)$ is the phase matching function, defined as

$$\mathrm{\Upsilon }\left({\omega }_1,{\omega }_2,{\omega }_3,{\omega }_4\right)={\int }_0^z{e}^{-i\left(k_{{\omega }_1}\left(x\right)-k_{{\omega }_2}\left(x\right)+k_{{\omega }_3}\left(x\right)-k_{{\omega }_4}\left(x\right)\right)}\mathit{dx}$$E32

Assuming the nonlinear Hamiltonian ${\widehat{H}}_1$ as a perturbative term of the Hamiltonian ${\widehat{H}}_0$ for which the continuous modes are noninteracting, and assuming the initial time of the interaction $t_0=-\infty$ and the final time $t_1=+\infty$, it is possible to relate the input field of the system to the output one introducing the asymptotic output field

$${\widehat{a}}_{\mathit{R\omega }}^{\mathit{\text{out}}}={\widehat{U}\widehat{a}}_{\mathit{R\omega }}{\widehat{U}}^{†},$$E33

where $\widehat{U}$ is the asymptotic unitary evolution operator (approximated to the first order in ${\widehat{H}}_1$)

$$\widehat{U}\equiv \widehat{U}\left(-\infty ,\infty \right)=e^{-\frac{i}{\mathrm{\hslash }}{\widehat{K}}_1}\text{where}{\widehat{K}}_1={\widehat{K}}_{\mathit{CPM}}+{\widehat{K}}_{\mathit{SQ}}+K_{\mathit{SPM}}$$E34

Working in the monochromatic degenerate pump limit $b\left({\mathrm{\Omega }}_{p^{\prime }}\right)=b\left({\mathrm{\Omega }}_p\right)\to b\left({\mathrm{\Omega }}_p\right)\delta \left({\mathrm{\Omega }}_p-{\omega }_p\right)\equiv b_p$ (where ${\omega }_p$ is the pump frequency), the propagators take the form

$${\widehat{K}}_{\mathit{CPM}}=-2\mathrm{\hslash }{\left|{\beta }_p\right|}^{2}z{\int }_0^{\infty }{k}_{{\omega }_s}{\widehat{a}}_{R{\omega }_s}^{†}{\widehat{a}}_{R{\omega }_s}d{\omega }_s$$E35

$${\widehat{K}}_{\mathit{SQ}}=-\frac{i\mathrm{\hslash }}{2}{\left|{\beta }_p\right|}^2{\int }_0^{\infty }\sqrt{k_{{\omega }_s}{k}_{{\omega }_i}}\left[\frac{1}{\mathrm{\Delta }k}\left(e^{i\mathrm{\Delta }\mathit{kz}}-1\right)\right]{\widehat{a}}_{R{\omega }_s}^{†}{\widehat{a}}_{R{\omega }_i}^{†}d{\omega }_s+\mathrm{H}.\mathrm{c}.$$E36

and

$$K_{\mathit{SPM}}=-\mathrm{\hslash }z{\left|{\beta }_p\right|}^2{k}_{{\omega }_p}{b}_p^{\ast }{b}_p$$E37

where $\beta \left({\omega }_p\right)\equiv {\beta }_p$ and $\mathrm{\Delta }k=2k_{{\omega }_p}-k_{{\omega }_s}-k_{{\omega }_i}$ are the chromatic dispersions.

Similar to the previous classical treatment, one can introduce the co-rotating framework by replacing the field operators with

$${\widehat{a}}_{R{\omega }_s}\left(z\right)={\widehat{\tilde{a}}}_{R{\omega }_s}\left(z\right)e^{2i{\left|{\beta }_p\right|}^2{k}_{{\omega }_s}z}\text{and}{b}_p\left(z\right)={\tilde{b}}_p\left(z\right)e^{i{\left|{\beta }_p\right|}^2{k}_{{\omega }_p}z}$$E38

In this framework, one can derive the following differential equation

$$\frac{\partial {\widehat{\tilde{a}}}_{R{\omega }_s}\left(z\right)}{\partial z}=\frac{i}{\mathrm{\hslash }}\left[\frac{d}{\mathit{dz}}{\widehat{K}}_1,{\widehat{\tilde{a}}}_{R{\omega }_s}\right]=i{\beta }_p^2\sqrt{k_{{\omega }_s}{k}_{{\omega }_i}}{e}^{i{\mathrm{\Psi }}_2\left({\omega }_s\right)z}{\widehat{\tilde{a}}}_{R{\omega }_i}^{†}$$E39

and

$$\frac{\partial {\tilde{b}}_p\left(z\right)}{\partial z}=-\left\{\frac{{\mathit{dK}}_{\mathit{SPM}}}{\mathit{dz}},{\tilde{b}}_p\right\}=0$$E40

where ${\mathrm{\Psi }}_2\left({\omega }_s\right)=\mathrm{\Delta }k+2{\left|{\beta }_p\right|}^2\left(k_{{\omega }_p}-k_{{\omega }_s}-k_{{\omega }_i}\right)$ is the total phase mismatch. These latter are formally identical to Eqs. (12) and (13), up to a frequency-dependent normalization of the wave-amplitudes, under the undepleted pump approximation. Reference [11] derives an exact solution for Eq. (39), being

$$\begin{array}{l}{\widehat{\tilde{a}}}_{R{\omega }_s}^{\mathit{\text{out}}}\left(z\right)=\left[\left(cosh\left(g_2\left({\omega }_s\right)z\right)-\frac{i{\mathrm{\Psi }}_2\left({\omega }_s\right)}{2g_2\left({\omega }_s\right)}sinh\left(g_2\left({\omega }_s\right)z\right)\right){\widehat{\tilde{a}}}_{R{\omega }_s}\right.\\ +i\frac{{\beta }_p^2\sqrt{k_{{\omega }_s}{k}_{{\omega }_i}}}{g_2\left({\omega }_s\right)}sinh\left(g_2\left({\omega }_s\right)z\right){\widehat{\tilde{a}}}_{R{\omega }_i}^{†}]e^{i\frac{{\mathrm{\Psi }}_2\left({\omega }_s\right)}{2}z}\end{array}$$E41

where

$$g_2\left({\omega }_s\right)=\sqrt{{\left|{\beta }_p^2\right|}^2{k}_{{\omega }_s}{k}_{{\omega }_i}-{\left(\frac{{\mathrm{\Psi }}_2\left({\omega }_s\right)}{2}\right)}^2}$$E42

If a state moves inside a TWJPA of length $z=\mathit{\text{aN}}$, the power gain will be $G\left({\omega }_s,\mathit{\text{aN}}\right)=⟨{\widehat{\tilde{a}}}_{R{\omega }_s}^{\mathit{\text{out}}}\left(\mathit{\text{aN}}\right)|{\widehat{\tilde{a}}}_{R{\omega }_s}^{{\mathit{\text{out}}}^{†}}\left(\mathit{\text{aN}}\right)⟩/⟨{\widehat{\tilde{a}}}_{R{\omega }_s}|{\widehat{\tilde{a}}}_{R{\omega }_s}^{†}⟩$.

3.3. Quantum Hamiltonian model based on discrete-mode operators

An alternative approach for the derivation of the quantum dynamics of a TWJPA is the one proposed in [32]. In this model, the quantum Hamiltonian for a 4WM parametric amplifier is expressed as the integral, along an arbitrary quantization length $l_q$, of the linear energy density stored in each element of the circuit. The energy stored per unit length in a Josephson junction can be derived from Eq. (5) by dividing each term by the elementary cell length $a$ and replacing ${\mathrm{\Delta \Phi }}_n$ with its continuous counterpart. Instead, the energy stored per unit length in a capacitance C can be alternatively expressed in terms of flux difference $\mathrm{\Delta \Phi }$ or stored charge Q as

$$U_C=\frac{1}{a}{\int }_{t_0}^t\mathit{VI}d{t}^{\prime }=\frac{1}{a}{\int }_{t_0}^t\frac{d\mathrm{\Delta \Phi }}{d{t}^{\prime }}C\frac{d}{d{t}^{\prime }}\left[\frac{d\mathrm{\Delta \Phi }}{d{t}^{\prime }}\right]d{t}^{\prime }=\frac{1}{2}\frac{C}{a}\mathrm{\Delta \Phi }\left(t\right)\frac{{\partial }^2\mathrm{\Delta \Phi }}{\partial t^2}$$E43

$$=\frac{1}{a}{\int }_{t_0}^t\frac{Q}{C}\frac{\mathit{dQ}}{d{t}^{\prime }}d{t}^{\prime }=\frac{1}{2a}\frac{1}{C}{Q}^2\left(t\right)$$E44

under the assumption that $\mathrm{\Delta \Phi }\left(t_0\right)=0$ and $Q\left(t_0\right)=0$.

Therefore, the quantum Hamiltonian of the system can be expressed, with an approximation to the first nonlinear order, as

$$\begin{array}{l}\widehat{H}={\int }_{l_q}\left[U_J+U_{C_J}+U_{C_0}\right]\mathit{dx}\approx \\ \approx {\int }_{l_q}\left[\left(\frac{1}{2{\mathit{aL}}_{J_0}}\mathrm{\Delta }\widehat{\mathrm{\Phi }}-\frac{1}{24{\mathit{aL}}_{J_0}}{\left(\frac{2\pi }{{\mathrm{\Phi }}_0}\right)}^2\mathrm{\Delta }{\widehat{\mathrm{\Phi }}}^3+\frac{1}{2}\frac{C_J}{a}\frac{{\partial }^2\mathrm{\Delta }\widehat{\mathrm{\Phi }}}{\partial t^2}\right)\mathrm{\Delta }\widehat{\mathrm{\Phi }}+\frac{1}{2a}\frac{1}{C_0}{\widehat{Q}}_{C_0}^2\right]\mathit{dx}\end{array}$$E45

where $\widehat{Q}$ and $\widehat{\mathrm{\Phi }}$ are quantum operators. The former can be expressed, as suggested in [34] and adapted for discrete-mode operators in [36], as

$${\widehat{Q}}_{C_0}=\sum _n\frac{C_0}{a}{\widehat{V}}_{C_0,n}=\sum _n\frac{C_0}{a}\sqrt{\frac{\mathrm{\hslash }{\omega }_{n}a}{2C_0{l}_q}}\left({\widehat{a}}_n{e}^{i\left(k_{n}x-{\omega }_{n}t\right)}+\mathrm{H}.\mathrm{c}.\right)$$E46

(here $k_{-n}=-k_n$ and ${\omega }_{-n}={\omega }_n$).

Before defining the flux operator, it is necessary to define an effective inductance $L_{\mathit{eff}}$ of the transmission line (modeled, as shown in Figure 2, as a parallel of the nonlinear Josephson inductance $L_J$ and the capacitance $C_J$):

$$\frac{1}{j{\omega }_n{L}_{\mathit{eff}}}=\frac{1}{j{\omega }_n{L}_J}+j{\omega }_n{C}_J\text{hence}{L}_{\mathit{eff}}=\frac{L_J}{1-{\omega }^2{L}_J{C}_J}\equiv L_J{\mathrm{\Lambda }}_n$$E47

Using the telegrapher’s equation [37], the discrete-mode current operator, under slowly varying amplitude approximation ($\partial {\widehat{a}}_n/\partial x\approx 0$), can be derived from the discrete-mode voltage operator as:

$${\widehat{I}}_{L_{\mathit{eff}}}=\sum _{n}sgn\left(n\right)\sqrt{\frac{\mathrm{\hslash }{\omega }_{n}a}{2L_{\mathit{eff}}{l}_q}}\left({\widehat{a}}_n{e}^{i\left(k_{n}x-{\omega }_{n}t\right)}+\mathrm{H}.\mathrm{c}.\right)$$E48

Therefore, the flux operator can be expressed as

$$\mathrm{\Delta }\widehat{\mathrm{\Phi }}=\frac{L_J}{a}{\widehat{I}}_{L_{\mathit{eff}}}\text{where}{L}_J\left(\mathrm{\Delta \Phi }\right)=\frac{\mathrm{\Delta \Phi }}{I_J}=\frac{\frac{2\pi }{{\mathrm{\Phi }}_0}\mathrm{\Delta \Phi }}{I_c\frac{2\pi }{{\mathrm{\Phi }}_0}sin\left(\frac{2\pi }{{\mathrm{\Phi }}_0}\mathrm{\Delta \Phi }\right)}=L_{J_0}\frac{\frac{2\pi }{{\mathrm{\Phi }}_0}\mathrm{\Delta \Phi }}{sin\left(\frac{2\pi }{{\mathrm{\Phi }}_0}\mathrm{\Delta \Phi }\right)}$$E49

The recursive relation deriving from Eq. (49) can be solved iteratively. Exploiting a power series expansion of the sine function and considering just the first order of interaction, it results that

$$\mathrm{\Delta }\widehat{\mathrm{\Phi }}=\sum _n[\left[1+\frac{{\mathrm{\Lambda }}_n}{12}{\left(\frac{2\pi }{{\mathrm{\Phi }}_0}\mathrm{\Delta }{\widehat{\mathrm{\Phi }}}^{\left(0\right)}\right)}^2+O\left[{\left(\frac{2\pi }{{\mathrm{\Phi }}_0}\mathrm{\Delta }{\widehat{\mathrm{\Phi }}}^{\left(0\right)}\right)}^4\right]\right]\mathrm{\Delta }{\widehat{\mathrm{\Phi }}}_n^{\left(0\right)}$$E50

where $\mathrm{\Delta }{\widehat{\mathrm{\Phi }}}^{\left(0\right)}$ is the zero-order approximation of the flux quantum operator

$$\mathrm{\Delta }{\widehat{\mathrm{\Phi }}}^{\left(0\right)}=\sum _n\mathrm{\Delta }{\widehat{\mathrm{\Phi }}}_n^{\left(0\right)}=\sum _n\frac{k_{n}a}{{\omega }_n}\sqrt{\frac{\mathrm{\hslash }{\omega }_{n}a}{2C_0{l}_q}}\left({\widehat{a}}_n{e}^{i\left(k_{n}x-{\omega }_{n}t\right)}+\mathrm{H}.\mathrm{c}.\right)$$E51

Substituting Eqs. (46) and (50) in Eq. (45), and limiting the expression to the first nonlinear order, the Hamiltonian for a 4WM amplifiers turns up

$$\begin{array}{l}\widehat{H}=\sum _n\mathrm{\hslash }{\omega }_n\left({\widehat{a}}_n^{†}{\widehat{a}}_n+\frac{1}{2}\right)+\sum _{n,m,l,k}\frac{-i{\mathrm{\hslash }}^{2}a}{96L_{J_0}{I}_c^2{l}_q^2\mathrm{\Delta }{k}_{\mathit{\text{nmlk}}}}{e}^{-i\mathrm{\Delta }{\omega }_{\mathit{\text{nmlk}}}t}\left(e^{i\mathrm{\Delta }{k}_{\mathit{\text{nmlk}}}{l}_q}-1\right)\\ \times \{\left(1-4L_{J_0}{\mathrm{\Lambda }}_n{C}_J{\omega }_k^2\right){\left(\widehat{\tilde{a}}+\mathrm{H}.\mathrm{c}.\right)}_{n\times m\times l\times k}+4L_{J_0}{\mathrm{\Lambda }}_n{C}_J[2{\left(\omega \left(-i\widehat{\tilde{a}}+\mathrm{H}.\mathrm{c}.\right)\right)}_{n\times m}\\ \times {\left(\widehat{\tilde{a}}+\mathrm{H}.\mathrm{c}.\right)}_{l\times k}+{\left(\widehat{\tilde{a}}+\mathrm{H}.\mathrm{c}\right)}_{n\times m}{\left(\omega \left(-i\widehat{\tilde{a}}+\mathrm{H}.\mathrm{c}.\right)\right)}_{l\times k}\left]\right\}\end{array}$$E52

where $\widehat{\tilde{a}}\equiv sgn\left(n\right)\sqrt{{\mathrm{\Lambda }}_n{\omega }_n}{\widehat{a}}_n$, $\mathrm{\Delta }{k}_{\mathit{\text{nmlk}}}\equiv \pm k_n\pm k_m\pm k_l\pm k_k$, $\mathrm{\Delta }{\omega }_{\mathit{\text{nmlk}}}\equiv \pm {\omega }_n\pm {\omega }_m\pm {\omega }_l\pm {\omega }_k$ (a ± sign refers to a corresponding annihilation (creation) operator) and the subscript $i\times j$ indicates a multiplication (i.e., ${\left(\mathrm{\Lambda }\omega \right)}_{i\times j}={\mathrm{\Lambda }}_i{\omega }_i{\mathrm{\Lambda }}_j{\omega }_j$).

Neglecting the constant zero-point energy and assuming a strong degenerate classical pump (as shown in [36])

$${\widehat{a}}_p\to -i\sqrt{\frac{{\omega }_p{C}_0{l}_q}{2\mathrm{\hslash }a}}{A}_p$$E53

it is possible to approximate the Hamiltonian in Eq. (52) to the second order in ${\widehat{a}}_{s,i}^{\left(†\right)}$ as

$${\widehat{H}}^{\left(\mathit{CP}\right)}\approx \sum _{n=s,i}\mathrm{\hslash }\left({\omega }_n+{\xi }_n{\left|A_p\right|}^2\right){\widehat{a}}_n^{†}{\widehat{a}}_n-\mathrm{\hslash }\left(\chi A_p^2{\widehat{a}}_s^{†}{\widehat{a}}_i^{†}+\mathrm{H}.\mathrm{c}.\right)$$E54

where

$${\xi }_n=\frac{k_p^2{a}^2{\mathrm{\Lambda }}_n{\omega }_n}{32I_c^2{L}_{J_0}^2}\left(4-3{\delta }_{\mathit{pn}}\right)\left(1+\frac{2}{3}\left(\frac{{\mathrm{\Lambda }}_p}{{\mathrm{\Lambda }}_n}+\frac{{\mathrm{\Lambda }}_n}{{\mathrm{\Lambda }}_p}-2\right)\right)$$E55

represents the quantum self-phase modulation (when $n=p$) and the quantum cross-phase modulation (when $n=s,i$), whereas the coupling constant $\chi$ is defined as

$$\begin{array}{l}\chi =\frac{k_p^2{a}^2\sqrt{{\mathrm{\Lambda }}_s{\omega }_s{\mathrm{\Lambda }}_i{\omega }_i}}{16I_c^2{L}_{J_0}^2}\left(1+\frac{L_{J_0}{C}_J}{6}\left[{\omega }_p{\omega }_s\left(-2{\mathrm{\Lambda }}_p+5{\mathrm{\Lambda }}_s-3{\mathrm{\Lambda }}_i\right)\right.\right.\\ +{\omega }_p{\omega }_i\left(-2{\mathrm{\Lambda }}_p-3{\mathrm{\Lambda }}_s+5{\mathrm{\Lambda }}_i\right)+{\omega }_s{\omega }_i\left(4{\mathrm{\Lambda }}_p-2{\mathrm{\Lambda }}_s-2{\mathrm{\Lambda }}_i\right)\left]\right)\end{array}$$E56

Starting from the Hamiltonian ${\widehat{H}}^{\left(\mathit{CP}\right)}$, it is possible to calculate the Heisenberg equation of motion for the classical pump amplitude and for the quantum operators ${\widehat{a}}_s$ and ${\widehat{a}}_i$, obtaining the coupled mode equations:

$$\frac{\partial A_p}{\partial t}=-i\left({\omega }_p+2{\xi }_p{\left|A_p\right|}^2\right)A_p+2i{\chi }^{\ast }{A}_p^{\ast }{\widehat{a}}_s{\widehat{a}}_i$$E57

$$\frac{\partial {\widehat{a}}_{s\left(i\right)}}{\partial t}=-i\left({\omega }_{s\left(i\right)}+{\xi }_{s\left(i\right)}{\left|A_p\right|}^2\right){\widehat{a}}_{s\left(i\right)}+\mathit{i\chi }{A}_p^2{\widehat{a}}_{i\left(s\right)}^{†}$$E58

In [32], the hypothesis under which the classical coupled mode equations Eqs. (12) and (13) can be obtained from Eqs. (57) and (58) is described in detail.

Moving to a co-rotating frame (${\widehat{a}}_{s\left(i\right)}\to {\widehat{a}}_{s\left(i\right)}{e}^{i{\xi }_{s\left(i\right)}{\left|A_{p_0}\right|}^{2}z}$), the Hamiltonian (Eq. (54)) can be expressed as

$${\widehat{H}}_{\mathit{rot}}^{\mathit{CP}}=-\mathrm{\hslash }\left(\chi {\left|A_p\right|}^2{\widehat{a}}_s^{†}{\widehat{a}}_i^{†}{e}^{-i{\mathrm{\Psi }}_3^{\text{'}}t}+\mathrm{H}.\mathrm{c}.\right)$$E59

where ${\mathrm{\Psi }}_3^{\prime }=\left(4{\xi }_p-{\xi }_s-{\xi }_i\right){\left|A_{p_0}\right|}^2$. In this frame, introducing the undepleted pump assumption (Eq. (58)) turns into

$$\frac{\partial {\widehat{a}}_{s\left(i\right)}}{\partial t}=\mathit{i\chi }{\left|A_{p_0}\right|}^2{\widehat{a}}_{i\left(s\right)}^{†}{e}^{-i{\mathrm{\Psi }}_3^{\text{'}}t}$$E60

whose solutions are

$${\widehat{a}}_{s\left(i\right)}\left(t\right)=\left[{\widehat{a}}_{s{\left(i\right)}_0}\left(cosh\left(g_3^{\text{'}}t\right)+\frac{i{\mathrm{\Psi }}_3^{\text{'}}}{2g_3^{\text{'}}}sinh\left(g_3^{\text{'}}t\right)\right)+\frac{\mathit{i\chi }{\left|A_{p_0}\right|}^2}{g_3^{\text{'}}}{\widehat{a}}_{i{\left(s\right)}_0}^{†}sinh\left(g_3^{\text{'}}t\right)\right]$$E61

where the exponential complex gain factor is defined as

$$g_3^{\prime }=\sqrt{{\left|\chi \right|}^2{\left|A_{p_0}\right|}^4-{\left(\frac{{\mathrm{\Psi }}_3^{\prime }}{2}\right)}^2}$$E62

If a state spends a time $t$ in the amplifier, the gain can be expressed as $G_s^Q\left(t\right)=⟨{\widehat{a}}_s^{†}\left(t\right){\widehat{a}}_s\left(t\right)⟩/⟨{\widehat{a}}_{s_0}^{†}{\widehat{a}}_{s_0}⟩$.

To make the results of this last treatment, in which the operators are expressed as a function of the time, comparable with the previous ones, in which the operators are expressed as a function of the space coordinate, we need to take into account the phase velocity of the tones. It turns out that:

$${\mathrm{\Psi }}_3=\mathrm{\Delta }k+\left(4{\xi }_p\frac{{\omega }_p}{\mid k_p\mid }-{\xi }_s\frac{{\omega }_s}{\mid k_s\mid }-{\xi }_i\frac{{\omega }_i}{\mid k_i\mid }\right)\text{and}{g}_3=\sqrt{{\left|\chi \right|}^2{\left|A_{p_0}\right|}^4{\left(\frac{{\omega }_p}{\mid k_p\mid }\right)}^2-{\left(\frac{{\mathrm{\Psi }}_3}{2}\right)}^2}$$E63

where $\mathrm{\Delta }k=2k_p-k_s-k_i$ is the chromatic dispersion.

3.4. Models comparison

The three exponential complex gain factors ($g_i$) and the three total phase mismatches (${\mathrm{\Psi }}_i$) derived in these models are analytically different but numerically similar, as shown in Figure 3 (where the two insets report the differences between the quantum predictions and the classical ones).

Figure 3.

Comparison of total phase mismatches (${\mathit{\Psi }}_{\mathit{i}}$) and exponential complex gain factors (${\mathit{g}}_{\mathit{i}}$) predicted by the three treatments presented in Section 3. For these calculations, typical constructive parameters have been assumed: $\mathit{a}=\mathit{50}\mathit{\mu m}$, ${\mathit{I}}_{\mathit{c}}=\mathit{5}\mathit{\mu A}$, ${\mathit{C}}_{\mathit{J}}=\mathit{300}\mathit{fF}$, and ${\mathit{C}}_{\mathit{0}}=\mathit{35}\mathit{fF}$, in such a way that the characteristic impedance is $\mathit{Z}\approx \mathit{50}\mathit{\Omega }$. Furthermore ${\mathit{\omega }}_{\mathit{p}}/\mathit{2}\mathit{\pi }=\mathit{6}\mathit{GHz}$ and ${\mathit{I}}_{\mathit{p}}=\mathit{Ic}/\mathit{2}$. In the insets, the differences between the quantum predictions and the classical ones are presented.

It is important to observe that in Eqs. (17), (41), and (61), the last term is always equal to zero in the case of a zero initial idler amplitude. In such a case, and under the hypothesis of a perfect phase matching (${\mathrm{\Psi }}_i=0$), $g_i$ is real and the amplification gain increases exponentially with the line length; whereas, in the case of a nonzero phase mismatch, $g_i$ is imaginary and the gain increases quadratically [11].

4. Conclusions

In the present chapter, we have presented the state-of-the-art of the experimental evidences in the field of Josephson junctions-based traveling-wave metamaterials through a historical review in Section 2. Moreover, in Section 3, we have reported three different theoretical approaches for the prediction of a TJWPA dynamics, in the particular case of a 4WM process. Assuming similar simplifying hypothesis, like the use of a classical undepleted degenerate pump, the presence of slowly varying fields along the transmission line and approximating the nonlinearities of the system up to the first order, a similar expression for the signal amplitude (or field annihilation/creation operators in the case of quantum theories) expressed in a co-rotating frame, is derived in the three treatments. Although the results of the quantum theories are similar to the classical ones, the description of the system dynamics with a quantum theory grants the possibility to evaluate photon-number distributions, squeezing effects and averages, standard deviations or higher-order moments of the measurements operators, taking into account the commutation relations between operators explicitly. For instance, detailed calculations of the output state of a TWJPA in the case of a single-photon input state and in the case of a coherent input state are presented in [32].

Acknowledgments

The author would like to thank Luca Callegaro for the stimulating discussion. This work was partially funded by the Joint Research Project PARAWAVE of the European Metrology Programme for Innovation and Research (EMPIR). This project has received funding from the EMPIR program co-financed by the Participating States and from the European Unions’ Horizon 2020 research and innovation program.