A quantum model for rf-SQUIDs based metamaterials enabling 3WM and 4WM Travelling Wave Parametric Amplification
Angelo Greco, Luca Fasolo, Alice Meda, Luca Callegaro, Emanuele Enrico
AA quantum model for rf-SQUIDs based metamaterials enabling 3WM and 4WMTravelling Wave Parametric Amplification
Angelo Greco and Luca Fasolo
INRiM, Istituto Nazionale di Ricerca Metrologica, Strada delle Cacce 91, 10135 Torino, Italy andDepartment of Electronics and Telecommunications,PoliTo, Corso Castelfidardo 39, 10129 Torino, Italy
Alice Meda, Luca Callegaro, and Emanuele Enrico ∗ INRiM, Istituto Nazionale di Ricerca Metrologica, Strada delle Cacce 91, 10135 Torino, Italy (Dated: September 22, 2020)A quantum model for Josephson-based metamaterials working in the Three-Wave Mixing (3WM)and Four-Wave Mixing (4WM) regimes at the single photon level is presented. The transmissionline taken into account, namely Traveling Wave Josephson Parametric Amplifier (TWJPA), is abipole composed by a chain of rf-SQUIDs which can be biased by a DC current or a magnetic fieldin order to activate the 3WM or 4WM nonlinearities. The model exploits a Hamiltonian approachto analytically work out the time evolution both in the Heisenberg and interaction pictures. Theformer returns the analytic form of the gain of the amplifier, while the latter allows to recover theprobability distributions vs. time of the photonic populations, for multimodal Fock and coherentinput states. The dependence of the metamaterial’s nonlinearities is presented in terms of circuitparameters in a lumped model framework while evaluating the experimental conditions effects onthe model validity.
I. INTRODUCTION
Superconducting amplifiers are nowadays widely used forthe detection of single photons in several ranges of the elec-tromagnetic spectrum. From microwaves to X-rays thesedevices have shown unrivalled performances for what con-cerns quantum efficiency, resolving power and added noise,compared to their solid state counterparts [1–6]. The pe-culiar characteristics of superconducting materials allow toengineer highly performing resonators and cavities, char-acterised by quality factor of the order of ≈ [7–10].Indeed, resonator-based superconducting amplifiers show aquite high gain, in the range of
20 dB [9], however, they aresubjected to a limited bandwidth, fact that makes them un-suitable for the multiplexing required in complex systems.Traveling Waves Josephson Parametric Amplifiers (TWJ-PAs) and Kinetic Inductance Traveling Wave Amplifiers(KITs) promise to be appropriate devices for this aim inthe microwave regime, showing in principle valuable multi-plexing capabilities due to their wide bandwidth [11]. In-deed, it has been shown how the Four-Wave Mixing (4WM)induced in all the Kerr-like media allows amplifying verytiny signals over a several
GHz bandwidth with a mini-mum amount of added noise [11–14]. Nevertheless, recentpapers show that, enabling the Three-Wave Mixing (3WM)interaction, through the introduction of a quadratic non-linearity in the medium, could provide several benefits andexperimental simplifications for what concerns feasibilityand integration capabilities. In particular, a three-wave ∗ Corresponding author: [email protected] mixer generally requires a lower input pump power, aneasier output filtering and shows a higher dynamic range[15–17]. These distinctive characteristics make TWJPAsworking in 3WM excellent candidates for the readout ofquantum limited detectors (e.g., rf-SETs, rf-SQUIDs), bypreserving the quantum properties of their outputs [18–20].Moreover, a three-wave mixer can be a promising candidatefor the generation of heralded photons pairs, since it natu-rally enables Parametric Down Conversion (PDC) [21].In this framework we develop a quantum model, based onprevious theoretical works [22, 23], for a recently proposedTWJPA concept [16] covering both the 3WM and 4WMregimes. Previous classical descriptions in terms of elec-tromagnetic waves [16, 17] were limited to the high powerrange, completely neglecting any description of the light-matter interaction at the single photon level. Our theoryexploits circuit-QED techniques to model a TWJPA madeup of a chain of rf-SQUIDs capacitively shunted to ground.The proposed layout can be biased by a DC current or anexternally applied magnetic field in order to activate 3WMor 4WM of the microwave traveling modes. The quantumdescription allows to analytically treat important figures ofmerit of the amplifier like gain, squeezing and time evolu-tion of arbitrary quantum states at the single photon level.The main results of the paper are reported in Section II.In particular, Subsection II A reports the Hamiltonian infirst quantization formalism, based on the circuit modelof a nonlinear lossless transmission line. Then, Subsec-tion II B is dedicated to develop the theory through theoccupation number formalism, and a 3WM/4WM Hamil-tonian is found. A selection of modes follows in SubsectionII C, leading to model the 3WM/4WM quantum mechani- a r X i v : . [ c ond - m a t . s up r- c on ] S e p cal phenomena in the Heisenberg picture. Solving the dy-namics of the system (i.e., Langevin equations) allows toanalytically calculate the gain and squeezing capabilities ofthe amplifier. In Subsection II D the time evolution of Fockand coherent input states due to nonlinear interactions isanalytically treated and on these basis various examplesof photon statistics in the Fock space are calculated. TheDiscussion section III is focused on both the circuit pa-rameter constraints bounded to physical and experimentalrequirements (Subsections III A) and on a clear outline ofthe limits of validity of the model (Subsection III C andIII B). II. RESULTSA. rf-SQUIDs array embedded in a transmission line
The TWJPA recently proposed [16] and theoreticallyquantum mechanically treated in this paper can be mod-eled as an array of rf-SQUIDs embedded in a superconduct-ing transmission line. In the following, the Hamiltonian ofthe system will be derived as a function of its circuit param-eters. As represented in Figure 1, each elementary cell iscomposed by a superconducting loop containing a Joseph-son junction (with its associated capacitance C J and induc-tance L J ) and a geometrical inductance L g . Furthermore,each loop is coupled to ground through a capacitor C g . Thesystem taken into account is non-dissipative and, in sakeof simplicity, all the elementary cells are considered equal. The length of the elementary cell along the z -direction (i.e.,the propagating direction of the modes) is defined as a .In presence of an electromagnetic field, each of these cellsstores a certain amount of energy that can be expressedas a function of the conjugate coordinates ˆΦ and ˆ Q , thegeneralized magnetic flux and charge at a certain node, re-spectively. The total amount of energy can be computed asthe sum of the energy stored in each of its components (seeAppendix A). Moreover, being the system under analysisa repetition of identical elementary units, the total energystored in the whole medium can be expressed as the sumof the energy stored in each cell.Under the assumption that the differences between the ˆΦ (and ˆ Q ) of a couple of consecutive nodes are small enough,these quantities can be considered as continuous functionsof time and space (i.e., ˆΦ( z, t ) and ˆ Q ( z, t ) ). This consid-eration can also be applied to the definition of the totalenergy stored in the medium, leading to an expression interms of an integral of a linear density Hamiltonian ˆ H : ˆ H = ˆ H L g + ˆ H L J + ˆ H C J + ˆ H C g (1)where in the right-hand side of equation (1) one can rec-ognize respectively the energy density associated to the ge-ometrical inductance L g , the Josephson inductance L J , theJosephson capacitance C J and the ground capacitance C g .Hence, the Hamiltonian of the whole system is the spatialintegral across the z-axis and along an interaction length l q = aN [24], of the linear Hamiltonian density in equation(1): ˆ H = (cid:90) l q ˆ H d z = 12 a (cid:90) l q d z I c ϕ (cid:32) − cos (cid:32) ∆ ˆΦ ϕ (cid:33)(cid:33) + 1 L g ∆ ˆΦ + C J (cid:32) ∂ ∆ ˆΦ ∂t (cid:33) + ˆ Q C g C g (2)having identified ∆ ˆΦ as the magnetic flux differencefunction across the cells, N as the number of unit cellscomposing the transmission line, I c as the critical currentof the Josephson junction, ϕ = (cid:126) / e as the reduced fluxquantum, being (cid:126) the reduced Plank constant and e theelementary charge.The presence of an external magnetic field or a DC cur-rent through the line induces a constant component in theflux difference across a cell. This means that ∆ ˆΦ can beconsidered as the sum of two components, a constant one ∆Φ DC and a time-dependent one δ ˆΦ∆ ˆΦ = ∆Φ DC + δ ˆΦ (3) B. Second quantization framework
Here the Hamiltonian will be expressed in terms of ladderoperators. In this view, the voltage drop on the groundcapacitors C g can be expressed using a mode decompositionassuming that sinusoidal waves are passing through the line[25] ˆ V C g ( z, t ) = (cid:88) n (cid:115) (cid:126) ω n C g N (cid:16) ˆ a n e i ( k n z − ω n t ) + H . c . (cid:17) (4)where ω n and k n are the angular frequency and wavenum-ber of the n -th mode while ˆ a is its annihilation operator.Positive indexes denote progressive waves ( k n > and ω n > ), while negative indexes denote regressive waves( k − n = − k n < and ω − n = ω n ).The link between the voltage drop and the current pass- a C J C J C g C g C g L J L J L g z ΔΦ = ΔΦ DC + δΦ V Cg Ĩ = I + I DC B L g C J L J L g C g Figure 1. Electrical equivalent of a repetition of three elemen-tary cells (periodicity a ) of the rf-SQUIDs based TWJPA. Eachcell consists of a superconducting loop containing a geometricalinductance L g , a Josephson junction, with an associated capac-itance C J and inductance L J , and a ground capacitor C g . Theseries can be biased both through an external DC magnetic field B and a flowing current I DC . ∆ ˆΦ is the magnetic flux differenceacross the nodes of a cell, while ˆ V C g is the voltage drop acrossthe ground capacitor. ing through a cell is straightforwardly found recalling theTelegrapher’s equation, that exploits the inductance of thecell for the n -th mode ˆ L n ∂ ˆ V n ∂z = − ˆ L n a ∂ ˆ I n ∂t (5) ˆ L n can be calculated as the parallel between the effectiveinductance L eff ,n (composed by the Josephson capacitance C J and the geometrical inductance L g , see Appendix B)and the nonlinear Josephson inductance ˆ L J . Exploiting theconstitutive relation for a generic inductor it can be writtenthat ∆ ˆΦ = ˆ L ˆ I . Hence, using the flux-current relation ofa Josephson junction, ˆ I J = I c sin (cid:16) ∆ ˆΦ /ϕ (cid:17) , the nonlinearJosephson inductance ˆ L J can be simply expressed as ˆ L J = ∆ ˆΦˆ I J = ϕ I c ∆ ˆΦ /ϕ sin (cid:16) ∆ ˆΦ /ϕ (cid:17) ≡ L J ∆ ˆΦ /ϕ sin (cid:16) ∆ ˆΦ /ϕ (cid:17) (6)with L J = ϕ /I c . It follows that the cell inductance L n can be written as ˆ L n = Λ n L g n L g L J0 sin ( ∆ˆΦ /ϕ )( ∆ˆΦ /ϕ ) (7)where the dispersion coefficient of the n -th node Λ n =1 / (1 − ω n L g C J ) (Appendix B) has been defined. The time-dependent component of equation (3) can be found exploit-ing the mode decomposition for the AC current through thecell ˆ I n and the inductance ˆ L n for the corresponding modeas δ ˆΦ = (cid:88) n ˆ L n ˆ I n (8) It follows that (see Appendix C) δ ˆΦ = (cid:88) n n L g L J sin (cid:16) ∆Φ DC + δ ˆΦ ϕ (cid:17) ∆Φ DC + δ ˆΦ ϕ − δ ˆΦ (0) n (9)where the zero order AC flux component of the n -thmode δ ˆΦ (0) n has been defined. Equation (9) is a recursiverelation for the nonlinear flux operator δ ˆΦ , which can bestraightforwardly solved at zero order by the substitution δ ˆΦ (cid:55)→ δ ˆΦ (0) in the right-hand side.In order to find an analytical solution one can perform theTaylor expansion of the square root of equation (9) and theJosephson energy into (3) for δ ˆΦ (0) (cid:28) ϕ . The maximumorder of expansion was chosen in order to take into ac-count scattering events involving at most 4 photons. Thisprocedure provides a valid approximation for the nonlineartime-dependent flux operator δ ˆΦ that can be substitutedinto equation (2) in order to obtain the Hamiltonian of thesystem in terms of ladder operators. ˆ H = (cid:126) χ + (cid:88) n (cid:126) χ ( n )1 (cid:18) ˆ a † n ˆ a n + 12 (cid:19) ++ (cid:88) n,l,m (cid:126) χ ( n,l,m )3 (cid:8) ˆ a + ˆ a † (cid:9) n,l,m δ ∆ ω n,l,m , ++ (cid:88) n,l,m,s (cid:126) χ ( n,l,m,s )4 (cid:8) ˆ a + ˆ a † (cid:9) n,l,m,s δ ∆ ω n,l,m,s , (10)The subscripts of the braces in equation (10) stand for amultiplication of the form { ˆ a + ˆ a † } n,l,...,k = (ˆ a n + ˆ a † n )(ˆ a l +ˆ a † l ) · ... · (ˆ a k + ˆ a † k ) . The δ ∆ ω, Kronecker functions have therole to select the only scattering events that fulfill the en-ergy conservation between the three ( ∆ ω n,l,m = 0 ) or four( ∆ ω n,l,m,s = 0 ) modes taken into account .With this in hand the full Hamiltonian of the system isfound to be composed by a sum of four terms, the last twobeing interaction terms, where three or more modes giverise to 3WM or 4WM. χ ( n )1 and χ describe respectively the free field energy ofthe traveling modes and the magnetic energy stored intothe rf-SQUIDs due to the magnetic field or DC currentbias applied. Furthermore χ ( n,l,m )3 and χ ( n,l,m,s )4 are re-spectively the coupling parameters that characterize the3WM and 4WM, both strongly dependent on the circuit ∆ ω n,l,m = ± ω n ± ω l ± ω m , where the sign of each addend is de-fined by the combination of creation and annihilation operatorsthat precede this quantity (minus sign if related to a creationoperator, plus sign if related to an annihilation operator). Forinstance (ˆ a † n e − i ( k n z − ω n t ) )(ˆ a † l e i ( k l z − ω l t ) )(ˆ a † m e − i ( k m z − ω m t ) ) =ˆ a † n ˆ a l ˆ a † m e i (∆ k n,l,m z − ∆ ω n,l,m t ) , where ∆ k n,l,m = − k n + k l − k m and ∆ ω n,l,m = − ω n + ω l − ω m Parameter Value Description I c µ A Josephson critical current C g
14 fF
Ground capacitance L g
53 pH
Geometrical inductance C J
60 fF
Josephson capacitance a µ m Unit cell length N
900 Number of unit cells ω p π ·
12 GHz
Pump frequency ω s π · Signal frequency ω i π · ω j π ·
19 GHz ∆Φ DC , /φ . ∆Φ DC , /φ . parameters of the unit cell and on the frequency of themodes that populate the TWJPA. For the complete ex-pression of the coupling coefficients as a function of thelayout and experimental parameters see Appendix D. Thedistinctive characteristic of the layout under study is thestrong dependence of these coupling parameters with re-spect to the external bias conditions, opening the possibil-ity to properly select a working regime (3WM or 4WM).Each coupling parameter, defined by a set of indices (e.g. n , l , m and s ), quantifies the interaction strength of the re-spective modes. It is then clear that different combinationsof indices represent different effects that take place into theTWJPA which contribute to the output field. If we nowfocus on a particular working regime of the amplifier it canbe noted that if the TWJPA is biased so that the Kerr-likenonlinearity is suppressed, it is legit to consider the ampli-fier as a pure three-wave mixer, where the conservation ofenergy imposes the creation of the so-called Idler mode atfrequency ω p − ω , being ω p and ω the Pump and Signalfrequencies respectively. On the contrary, if the quadraticnonlinearity gets suppressed, the TWJPA becomes a purefour-wave mixer, hence, with the degenerate pump assump-tion valid since now on for the 4WM regime, the Idler willbe located at ω p − ω . C. Gain and Squeezing in 3WM and 4WM
In order to analytically treat the problem the numberof traveling modes that populate the TWJPA will be re-stricted to three, the input pump and signal frequenciesplus the idler frequency that changes depending on theactive nonlinearity. This assumption implies that mixed3WM/4WM states will not be taken into account, becausethe presence of a mixed state would require a four coupledmodes discussion, that goes beyond the scope of this paper.Furthermore, since now the 3WM regime will be consideredin a non-degenerate condition, that is ω (cid:54) = ω p / . This factallows not to consider the phase difference of the incoming - χ i [ M H z ] χ · - ( χ ( p ) - ω p )· - ( χ ( i ) - ω i )· - ( χ ( s ) - ω s )· - ( χ { p ,s,i } )· - χ { p,p,s,j } - - ΔΦ DC / Φ ξ nn [ M H z ] ξ pp ξ ss ξ jj ξ ps ξ pj ξ sj χ { p,p,s,j } / ℏ / ℏ Figure 2. Hamiltonian coupling parameters χ i characteriz-ing the Hamiltonian (10) versus the normalized DC flux bias( ∆Φ DC / Φ ). The Hamiltonian coefficient related to the con-stant flux bias ( χ ) has been scaled by a factor of while thenon-interacting-modes Hamiltonian coupling constants ( χ ( n )1 )have been shifted by the frequency of the corresponding pho-ton ( ω i ) and scaled by a factor . The indices in the super-scripts vary with the considered mode and can take the values p(pump), s (signal) and i (idler). The blue vertical lines indicatethe flux biases at which the amplifier works as a three-wavemixer (see Section II C for a detailed description), while thered vertical lines indicate the flux biases at which the amplifierworks as a four-wave mixer. The coupling parameters ξ n,n and ξ n,l refer to the self-phase and cross-phase modulation due tothe 4WM interaction. The circuit parameters used to performthe numerical evaluations and plots are summarized in table I. waves in the following because the non-degenerate para-metric amplification results to be phase-preserving [22].With this in hand the full Hamiltonian (10), can be re-duced to two different forms depending on the regime theamplifier is working in. Concerning the 3WM, the followingHamiltonian is obtained ˆ H = (cid:126) χ + (cid:88) n = { ω p ,ω,ω p − ω } (cid:126) χ ( n )1 (cid:18) ˆ a † n ˆ a n + 12 (cid:19) ++ (cid:126) χ { ω p ,ω,ω p − ω } (cid:16) ˆ a † ω p ˆ a ω ˆ a ω p − ω + ˆ a † ω ˆ a † ω p − ω ˆ a ω p (cid:17) (11)having introduced χ { ω p ,ω,ω p − ω } as the sum of all the pos-sible terms arising from index permutations of χ ( ω p ,ω,ω p − ω )3 neglecting permutations signs degeneracy. While the 4WMHamiltonian results to be ˆ H = (cid:126) χ + (cid:126) ξ + (cid:88) n = { ω p ,ω, ω p − ω } (cid:126) χ ( n )1 (cid:18) ˆ a † n ˆ a n + 12 (cid:19) ++ (cid:88) n = { ω p ,ω, ω p − ω } (cid:126) ξ n ˆ a † n ˆ a n ++ (cid:88) n,l = { ω p ,ω, ω p − ω } (cid:126) ξ n,l ˆ a † n ˆ a n ˆ a † l ˆ a l ++ (cid:126) χ { ω p ,ω p ,ω, ω p − ω } ·· (cid:16) ˆ a † ω p ˆ a † ω p ˆ a ω ˆ a ω p − ω + ˆ a † ω ˆ a † ω p − ω ˆ a ω p ˆ a ω p (cid:17) (12)where ξ is a small correction to the zero-point energy, ξ n is a small contribution to the free-field energy of themodes and ξ n,l is the coefficient describing the self- ( n = l )and cross-phase ( n (cid:54) = l ) modulation phenomena. All theselatter arise from the 4WM interaction in the medium. Like-wise the 3WM case, χ { ω p ,ω p ,ω, ω p − ω } is the sum of allthe possible terms that derive from index permutations of χ ( ω p ,ω p ,ω, ω p − ω )4 neglecting permutations signs degeneracy.Figure 2 shows the behaviour of the most significant cou-pling parameters as a function of ∆Φ DC . All these coef-ficients present an oscillating-like behaviour with multiplezeros for different values of ∆Φ DC while red and blue ver-tical lines represent particular bias values (working points)that select the 4WM or 3WM working regimes respectively.Since now on we recall the 3WM and 4WM regimes by re-ferring to the first blue and red vertical lines in Figure 2.For their numerical values see Table I.Once the Hamiltonian of the system is known, it is pos-sible to work out the dynamic of the observables. Ex-ploiting the Heisenberg picture of quantum mechanic, thetime evolution of the creation and annihilation opera-tors can be computed through the Heisenberg equation dˆ a H ( t ) / d t = ( i/ (cid:126) )[ ˆ H, ˆ a H ( t )] + ( ∂ ˆ a/∂t ) H (for the completecalculations see Appendix E). From here to the end of Sec-tion II C we will drop the H subscript.A system of coupled equations (3WM (E1)-(E3), 4WM(E4)-(E6)) comes out, that is not exactly solvable ana-lytically [26]. Indeed, one can proceed with the so-called undepleted pump approximation to analytically treat thesystem. Such an approximation requests the pump ampli-tude to be much higher than the signal and idler ones sothat its magnitude does not change significantly during theinteraction process. On the other hand, under the so-called classical pump approximation , the ladder operator describ-ing the pump mode can be treated as a classical amplitude (cid:115) (cid:126) ω p C g N ˆ a p (cid:55)→ A p (13)being A p the classical voltage amplitude of ˆ V C g (equa-tion (4)). The strong interplay between the traveling waves mani-fests itself in a system of coupled differential equations forthe annihilation operators describing the signal and idlermodes dˆ a ω d t = − i Υˆ a † ω (cid:48) e − i Ψ t (14a) dˆ a ω (cid:48) d t = − i Υˆ a † ω e − i Ψ t (14b)where the density phase mismatch Ψ has been defined(3WM - equation (E22), 4WM - equation (E15)). Themain structure of the system stays the same regardless ofthe kind of interaction that takes place into the TWJPA,indeed it is possible to define an interaction parameter Υ = Υ , that characterizes the working regime theamplifier is biased in. Υ = χ | A p , | (15a) Υ = χ | A p , | (15b) χ , are two bias tunable coefficients that incorporate in-formation about the strength of the quadratic or cubicnon-linearity into the device (for their definition refer toequations (E16) and (E23)). It has to be noticed that in Υ , the proportionality to the initial pump ampli-tude A p , reflects the nature of the scattering taken intoaccount, hence involving one (linear) or two (quadratic)pump photons.Under the undepleted pump assumption and working inthe co-rotating frame one can find the following analyticalsolution to equations (14a) and (14b) ˆ a ω ( t ) = (cid:34) ˆ a ω, (cid:32) cosh ( gt ) + i Ψ2 g sinh ( gt ) (cid:33) −− i Υ g (ˆ a ω (cid:48) , ) † sinh ( gt ) (cid:35) e − i (Ψ / t (16)being ˆ a ω, and (ˆ a ω (cid:48) , ) † the ladder operators at the initialinteraction time and with the complex gain factor g = (cid:115) Υ − (cid:18) Ψ2 (cid:19) (17)For the 3WM case, under experimentally reasonable pa-rameter (see Table I) a negligible total phase mismatch ap-proximation , hence the phase mismatch density times theinteraction time, can be considered in equation (16), so that Ψ t ≈ and the phase lag between the traveling modes canbe neglected. Moreover, in the undepleted pump approx-imation it can be shown that the gain variation given bythe phase mismatch density in (17) can be neglected since Υ (cid:29) Ψ , giving the much simpler relation g ≈ | Υ | (18) It is helpful to introduce a set of auxiliary functions thatincorporates the behaviour of the TWJPA and simplifiesthe notation u ( ω, t ) = cosh ( g ( ω ) t ) + i Ψ( ω )2 g ( ω ) sinh ( g ( ω ) t ) (19) v ( ω, t ) = − Υ g ( ω ) sinh ( g ( ω ) t ) (20) It is now possible to define the gain of a TWJPA as theratio among the average number of photons of frequency ω after a certain amount of time t spent into the mediumand the average number of photons at the same frequencyentering the transmission line G ( ω ) = (cid:104) ˆ a † ω ˆ a ω (cid:105)(cid:104) (ˆ a ω, ) † ˆ a ω, (cid:105) = | u | + | v | (cid:2) (cid:104) (ˆ a ω (cid:48) , ) † ˆ a ω (cid:48) , (cid:105) + 1 (cid:3) + iu ∗ v (cid:104) (ˆ a ω, ) † (ˆ a ω (cid:48) , ) † (cid:105) − iuv ∗ (cid:104) ˆ a ω (cid:48) , ˆ a ω, (cid:105)(cid:104) (ˆ a ω, ) † ˆ a ω, (cid:105) (21)Equation (21) is a general relation to estimate the num-ber of outgoing signal photons regardless of the natureof the incoming state (Fock, coherent, thermal, etc.).The classical limit of equation (21) can be found bysetting a large amount of input signal photons, hence (cid:104) (ˆ a ω, ) † ˆ a ω, (cid:105) (cid:29) and by setting the number of input idlerphotons to zero, so (cid:104) (ˆ a ω (cid:48) , ) † ˆ a ω (cid:48) , (cid:105) ≈ G ( ω ) ≈ | u | ≈ gt (22)that results to be in accordance with [16]. Moreover,it is possible to explicitly compute the gain for an in-coming Fock state | Ψ (cid:105) = | n ω , n ω (cid:48) (cid:105) or a coherent state | Ψ (cid:105) = | α ω , β ω (cid:48) (cid:105) starting from equation (21) G F , H = | u | + n ω (cid:48) + 1 n ω | v | (23) G C , H = | u | + ( | β ω (cid:48) | + 1) | v | + iu ∗ vα ∗ ω β ∗ ω (cid:48) − iuv ∗ α ω β ω (cid:48) | α ω | (24)Figure 3(a) shows several plots of equation (23) calcu-lated for Fock input states with different numbers of signalphotons and for the experimental parameters reported inTable I. It is evident that the quantum gain of (23) tends to(22) for a large number of input photons ( n ω > ). More-over, a closer look at the curves representing equation (22)with and without negligible phase mismatch shows thatthe approximation Ψ ≈ holds in all the bandwidth onlyfor the 3WM regime. The orange, red and blue curves areplotted for growing pump powers, hence the gain, whichshows a maximum at half the pump frequency, is stronglydependent on the pumping power.The correlation of the signal and idler photons results ina squeezed output field of the TWJPA. In order to model these correlations one can introduce quadratures as ˆ Y θ ( ω ) = i ( e iθ/ ˆ a † ω − e − iθ/ ˆ a ω ) (25)with their associated fluctuations ∆ ˆ Y θ ( ω ) = ˆ Y θ ( ω ) − (cid:104) ˆ Y θ ( ω ) (cid:105) (26)being θ the so-called squeezing angle.From the previous definitions one can compute (see Ap-pendix F) the relation between the squeezing spectrum S and the quadratures fluctuations as S ( ω ) = (cid:88) n (cid:104) ∆ ˆ Y θ ( ω )∆ ˆ Y θ ( ω n ) (cid:105) (27)For a vacuum input state, it can be shown that the prod-uct of the fluctuations of the two quadratures gives theminimum possible value allowed by the Heisenberg uncer-tainty principle, fingerprint of a quantum limited amplifi-cation [27]. From (27) the squeezing spectrum is then S ( ω ) = 1 + 2 | v ( ω, t ) | − | v ( ω, t ) | (cid:112) | v ( ω, t ) | + 1 (28)Figure 3(b) shows the squeezing spectrum of equation(28) plotted as a function of the signal frequency for dif-ferent input pump powers, calculated for a vacuum inputstate. The squeezing spectrum shows a minimum at halfthe pump frequency, just like the gain, that become morepronounced with a stronger pump power. D. Interaction of quantum states through 3WM or4WM
The time evolution of the state vectors can give impor-tant hints on the behaviour of the TWJPA in presence ofsingle photon level signals. Moving to the framework of | >| >| >| >| >| > ψ = ψ ≠ ψ = ψ ≠ G [ d B ] ( a ) - - - - - - ω s [ GHz ] S [ d B ] ( b ) Figure 3. (a) Gain G of the TWJPA under the undepletedpump approximation expressed by equation (22) and (23) inthe 4WM (light blue/orange) and 3WM (blue/red) regimes.Several curves with and without zero phase mismatch, respec-tively Ψ = 0 and Ψ (cid:54) = 0 , are presented. The lighter curvesare calculated considering a different number of input signalphotons for Fock states, while the purple and green indicatethe classical limit, expressed by equation (22) for I p /I c = 0 . ( P p = −
58 dBm ) and I p /I c = 1 . ( P p = −
55 dBm ) respectively.(b) Squeezing spectrum S (equation (28)) as a function of thesignal frequency calculated for a vacuum input state in the 4WMand 3WM regimes. Different colors express different pump cur-rents ( I p ) for which G and S are calculated: blue/light blue I p /I c = 0 . , red/orange I p /I c = 1 . . The working points in the3WM and 4WM cases are respectively ∆Φ DC , /φ = 0 . and ∆Φ DC , /φ = 0 . . the interaction picture it is possible to calculate the out-put photon statistics in the Fock base for any incomingstate. The time evolution of a quantum state | ψ ( t ) (cid:105) can beexpressed as | ψ ( t ) (cid:105) = e − i (cid:126) (cid:82) t ˆ H int d t (cid:48) | ψ (0) (cid:105) = e − i (cid:126) ˆ H int t | ψ (0) (cid:105) (29)where ˆ H int = ˆ H int,3WM(4WM) is the three(four)-wave mixing Hamiltonian written in the co-rotating frame, un-der the undepleted pump approximation: ˆ H int,3WM = (cid:126) χ | A p , | (cid:16) ˆ a s ˆ a i e i Ψ t + ˆ a † s ˆ a † i e − i Ψ t (cid:17) (30) ˆ H int,4WM = (cid:126) χ | A p , | (cid:16) ˆ a s ˆ a i e i Ψ t + ˆ a † s ˆ a † i e − i Ψ t (cid:17) (31)Under the negligible phase mismatch condition (i.e., Ψ t (cid:28) ) equation (29) becomes | ψ ( t ) (cid:105) = e iκ ( ˆ a s ˆ a i +ˆ a † s ˆ a † i ) | ψ (0) (cid:105) (32)where κ = − χ | A p , | t ( κ = − χ | A p , | t ) is the amplifi-cation factor for the 3WM (4WM) regime. Equation (32)can be written in a normal ordered form [28] as | ψ ( t ) (cid:105) = e i tanh ( κ )ˆ a † s ˆ a † i ·· e − ln [cosh ( κ )] ( a † s ˆ a s +ˆ a † i ˆ a i ) ·· e i tanh ( κ )ˆ a s ˆ a i | ψ (0) (cid:105) (33)In the following, the time evolution of two differentclasses of initial input states will be analyzed.
1. Fock States input
This subsection focus on the characteristics of the time-evolution of an initial Fock state | ψ F(0) (cid:105) = | N Sin (cid:105) s | N Iin (cid:105) i .Starting from equation (33) and considering their actionon the initial state, by means of a power expansion of eachexponential function, the expression of the quantum stateat a certain time t can be derived.Then, the expectation value of the signal photon numberoperator ˆ n s = ˆ a † s ˆ a s on the final state | ψ F ( t ) (cid:105) can be ex-pressed as (cid:104) ˆ n s (cid:105) ψ F ( t ) = (cid:104) ψ F ( t ) | ˆ n s | ψ F ( t ) (cid:105) = (cid:88) N Sfin P F (cid:0) N Sfin (cid:1) · N Sfin (34)where P F ( N Sfin ) is the probability to measure N Sfin sig-nal photons in the final state, and N Sin − min { N Sin , N
Iin } Iin } + 1 . This value reflects the number of possi-ble combinations that can occur between the initial signaland idler photons at the beginning of the interaction.Considering the case of an initial state | (cid:105) s | (cid:105) i (Fig. 4(b)),at t = 0 and with a certain probability, two couple ofsignal-idler photons may recombine to create a pair ofpump photons. This leads to the effective propagationand amplification of a single remaining signal photon.Yet, with a different given probability, just a singlecouple of signal-idler photons or none of them recombine,leading to the effective propagation and amplification of,respectively, two or three signal photons. In addition,considering Fig. 4(c) it can be noted that in the caseof N Sin = N Iin , despite the the fact that N Sin (cid:54) = 0 , theprobability to observe at the and of the amplifier N Sfin = 0 is significantly non-zero. This is in accordance with thefact that an effective propagation and amplification of thevacuum state can occur.Furthermore, exploiting the definition given in equation(34), and consistently with equation (23), a quantum defi- Log [ P C ( N finS )]- - - - - - N finS t / t T ( a ) N finS ( b ) N finS ( c ) Figure 5. Time evolution inside the medium, from its inputport ( t = 0 ) to the output port ( t = t T ) of the probability dis-tribution P C to find N Sfin signal photons for three different initialbimodal coherent states | α (cid:105) s | β (cid:105) i . (a) | (cid:105) s | (cid:105) i , (b) | (cid:105) s | (cid:105) i , and | (cid:105) s | (cid:105) i . The dashed purple lines represent the time evolutionof the expectation value (cid:104) ˆ n s (cid:105) nition of the gain of the amplifier can be given G F , I = (cid:104) ˆ n s (cid:105) ψ F ( t ) (cid:104) ˆ n s (cid:105) ψ F(0) = (cid:104) ˆ n s (cid:105) ψ F ( t ) N Sin (36) 2. Coherent States input Similarly to what has been performed in the case of aFock state input, the expectation value of the signal pho-ton number operator can be derived by (34) consideringan initial bimodal coherent state | ψ c (0) (cid:105) = | α (cid:105) s | β (cid:105) i andhaving the following probability distribution P C = ∞ (cid:88) m,n,n (cid:48) =0 ( − n − n (cid:48) [tanh ( κ )] n + n (cid:48) [cosh ( κ )] N Sfin + m − n (cid:48) ) ·· α N Sfin − n ( α ∗ ) N Sfin − n (cid:48) β m ( β ∗ ) m + n − n (cid:48) e [ | α | + | β | + i ( α ∗ β ∗ − αβ ) tanh ( κ )] ·· m ! ( N Sfin − n )! (cid:18) N Sfin n (cid:19)(cid:18) m + nn (cid:48) (cid:19) (37)The time evolution of the probability distribution P C is presented in figure (5) for three different initial bimodalcoherent states. In contrast with P F , this distribution is al-ways single-peaked over the whole range of the interactionand its maximum shifts in time starting from N Sfin = | α | .It can also be noticed that, for a fixed α , the distributionbecomes wider and wider with the increase of β .In this case again, the quantum gain of the amplifierresults to be consistent with the equation found in theHeisenberg picture (24) and can be expressed as G C , I = (cid:104) ˆ n s (cid:105) ψ c ( t ) (cid:104) ˆ n s (cid:105) ψ c (0) = (cid:104) ˆ n s (cid:105) ψ c ( t ) | α | (38) III. DISCUSSIONSA. Experimental constrains and impedance matching In order to couple the TWJPA with its electromagneticenvironment a particular (e.g. Z c = 50 Ω ) impedancematching is commonly required. This target can be reachedwith non-trivial additional on-chip components or by prop-erly tuning the cells parameters. This subsection is devotedto the search for a set of parameters that fulfill this aimby expressing all the other cell parameters as a function of I c . For simplicity and without lacking of generality, in thefollowing C J is supposed to be constant.In order to keep the induced magnetic flux function intothe rf-SQUID single valued, a design characterized by ascreening parameter β , given by β = 2 πL g I c φ < (39)is required. It is evident that a certain β set a hyperbolicrelation between L g and I c .For a generic mode n , an expression for C g having set I c (consequently L g ) and Z c can be inferred starting fromthe relation for the characteristic impedance of a losslesstransmission line ( Z n = (cid:113) L n /C n g ) C n g = L n Z n = 1 Z n Λ n L g n L g L J = L g − L g C J ω n Z n (cid:18) − L g C J ω n L g ϕ I c ∆Φ /ϕ /ϕ (cid:19) (40)It has to be noticed that the impedance matching canbe achieved just for a single mode since the character-istic impedance Z n of the line is frequency dependant.The matched mode can be engineered ad-hoc dependingon the experiment requirements. If a low power reflectionis mandatory, the matched mode should be the pump, in-stead, if no signal loss is preferred the signal mode shouldbe the matched one.Figure 6(a) reports several curves representing the trendsgiven by equations (39) and (40) plotted as functions of I c for different values of β and for a 50 Ω matching of a signalat . ( a ) L g [ p H ] C g [ f F ] β = β = β = ( b ) I c [µ A ] G F , H [ d B ] Figure 6. (a) The plot shows three sets of curves calculated fordifferent values of the screening parameter β representing thecell parameters for a 50 Ω matching of the signal mode at as a function of the critical current I c . The full curves refer tothe left axis and reports the geometrical inductance ( L g ) vs. I c .On the contrary, the dashed curves refer to the right axis andrepresents the ground capacitance ( C g ) vs. I c . (b) The figureshows the gain G F , H of the TWJPA in 3WM mode (full curves)and 4WM mode (dashed curves) as a function of the criticalcurrent I c for different values of the screening parameter β . Theground capacitance C g and the geometrical inductance L g arefixed for every value of I c by the 50 Ω matching condition (40)and by the value of β (39). Below the grey curves given by (42)the undepleted pump approximation holds. B. Validity of the undepleted pump approximation In order to obtain equations (14a) and (14b) the un-depleted pump approximation has been made. Physicallyspeaking this means that from there on the pump powershould be considered much higher than the signal and idlerones. This boundary directly translate into a common re-lation for the number of pump photons into the amplifier,that must be always at least 10 times higher than the othermodes. n p > · n s , i (41)Equation (41) is nothing but a condition on the modespowers that has been recast using the number of photons.Exploiting the definition of gain (21) and substituting the0operators with their expectation values, (41) becomes n s < n p Gn s , < n p G < n p n s , (42)the latter being a simple but powerful constraint on themaximum gain that the amplifier can express remainingwell described by the undepleted pump approximation.This condition is very general because it just depends onthe ratio between the number of pump photons and thenumber of incoming signal photons setting the limit atwhich the signal photonic population reaches the same or-der of magnitude then the pump one, hence the rate ofannihilation of pump photons due to the generation of sig-nal ones is no longer negligible. This limit is represented infigure 6(b) by the grey curves that cut the gain for a singlephoton Fock input state given by (23). For the screeningparameter β approaching unity the gain functions loose va-lidity for higher values of the critical current because thenon-linearity related to the induced flux into the rf-SQUIDgets stronger. For lower values of β the range of validitygets extended and the gain functions remain valid for lowercritical currents. C. Model validity due to Taylor expansions In subsection II B two different Taylor expansions wereperformed regarding the Josephson energy into (3) and thenonlinear flux operator (9). In both expansions the phaseswing δ Φ is considered to be small, reflecting the amplitudeof the AC current flowing into the transmission line. Forthis reason, it is necessary to point out the limit of valid-ity of the model in terms of phase swing, hence of current.Since the pump current is considered to be much higher ofthe signal and idler ones (undepleted and classical pumpapproximations) it is legit to consider all the current flow-ing into the TWJPA equal to the pump current I p . Thisfact means that, for the model to be valid, I p needs to besmaller than a certain threshold.An error function can be built for the Josephson energyand equation (9) as the difference between the real func-tion and its series expansion. The threshold is then chosenso that the error functions are always smaller than , forany I p used during the computations.Figure 7 shows the error functions calculated for equa-tion (9) (solid lines) and (3) (dashed lines) for ∆Φ DC valuescorresponding to the 3WM and 4WM working points. Thisconstraint fixes the maximum pump current allowed by themodel. It’s worth noting that the maximum allowed I p value is greater than the Josephson critical current. Thisregime is accessible since the single cell is composed by - - - - R e l . E rr . I p / I c Figure 7. The solid curves represent the error function of(9) while the dashed curves indicate the error function of theJosephson energy into (3). The error functions are calculatedas the relative difference between the function and its Taylorexpansion versus the normalized AC pump ( I p /I c ) calculatedin the 4WM and 3WM working points. The horizontal grayline indicates the threshold chosen as the reference error. the parallel of a geometric inductance and a Josephsonjunction, hence the total pump current can split unevenlybetween these two components, keeping the current flow-ing into the Josephson junction always below the criticalvalue. It’s worth mentioning that the resonances in the4WM curves arise from the I p values for which the Joseph-son inductance diverges. IV. CONCLUSIONS A quantum theory for the parametric amplification viaa chain of rf-SQUIDs embedded in a in a wave-guide hasbeen developed through a circuit-QED approach. A mixedlumped/distributed-element approach has been adopted inorder to define the Hamiltonian of the system, valid forboth 3WM and 4WM interactions. The dynamics of thesystem has been calculated first in the Heisenberg pic-ture where, through the solution of a system of QuantumLangevin Equations for the traveling modes, a closed formfor the evolution of the photonic populations, power gainand squeezing spectrum were found. Then, by means ofthe interaction picture, the time evolution of some repre-sentative input states (Fock and coherent states) has beencalculated, allowing to model the quantum dynamics ofphotonic amplification and recombination into the TWJPAin the few photons regime. The closed forms derived for thetransmission line parameters represents simple constraintson the unit cell circuit components, so that the model va-lidity has been linked to a clear experimental parameterspace.1 DATA AVAILABILITY All data generated or analysed during this study are in-cluded in this published article. Appendix A: Hamiltonian linear density of theelementary cell components In this appendix we start describing the general methodused to calculate the energy stored in a circuit element,then we derive the energy stored in each element whichconstitutes an elementary cell of a TWJPA.Defining I as the current flowing through a certain circuitelement and V as the voltage drop across it, the energystored in the electrical component at a certain time t canbe expressed as the time-integrated power P = V I : U ( t ) = (cid:90) tt P ( t (cid:48) )d t (cid:48) = (cid:90) tt I ( t (cid:48) ) · V ( t (cid:48) )d t (cid:48) (A1)The current flowing through a generic inductance L in-duces a magnetic flux ∆Φ( t ) = LI ( t ) , and can be relatedto the voltage drop across the element by the relation V ( t ) = L d I ( t )d t (A2)Hence one can express the energy stored in the geomet-rical inductance L g as U L g ( t ) = (cid:90) tt I L g ( t (cid:48) ) · V L g ( t (cid:48) )d t (cid:48) = (cid:90) tt I L g ( t (cid:48) ) · L g d I L g d t (cid:48) d t (cid:48) = L g I L g ( t ) = L g (cid:18) ∆Φ( t ) L g (cid:19) = (∆Φ( t )) L g (A3)having assumed I L g ( t ) = 0 .Exploiting the relation between magnetic flux differenceand voltage drop V ( t ) = d∆Φ( t )d t (A4)the energy stored in the nonlinear Josephson inductance L J can be expressed as U L J ( t ) = (cid:90) tt I L J ( t (cid:48) ) · V L J ( t (cid:48) )d t (cid:48) = (cid:90) tt I c sin (cid:18) ∆Φ( t (cid:48) ) ϕ (cid:19) · d∆Φ( t (cid:48) )d t (cid:48) d t (cid:48) = ϕ I c (cid:18) − cos (cid:18) ∆Φ( t ) ϕ (cid:19)(cid:19) (A5)having assumed ∆Φ( t ) = 0 .Exploiting the relation between the current flowing througha capacitance C and the voltage drop V across its terminals I ( t ) = C d V ( t )d t (A6) the energy stored in the ground capacitance C g can beexpressed as U C g ( t ) = (cid:90) tt I C g ( t (cid:48) ) · V C g ( t (cid:48) )d t (cid:48) = (cid:90) tt C g d V C g ( t (cid:48) )d t (cid:48) · V C g ( t (cid:48) )d t (cid:48) = C g V C g ( t ) = 12 C g Q C g (A7)having assumed V C g ( t ) = 0 .Lastly, exploiting relations (A4) and (A6), the energystored in the capacitance associated to the Josephson junc-tion can be expressed as U C J ( t ) = (cid:90) tt I C J ( t (cid:48) ) · V C J ( t (cid:48) )d t (cid:48) = C J (cid:90) tt dd t (cid:48) (cid:20) d∆Φ( t (cid:48) )d t (cid:48) (cid:21) · d∆Φ( t (cid:48) )d t (cid:48) d t (cid:48) = C J (cid:18) d∆Φ( t )d t (cid:19) (A8)having assumed ∆Φ( t ) = 0 .Using a standard procedure [25], one can derive theHamiltonian operator that describes an electrical circuitstarting from the definition of the energy stored in each ofits components and transforming the physical observablesinto the corresponding operators. Furthermore, being a asthe unit cell length, one can express the linear density ofHamiltonian associated to each component of the circuitrepresented in Figure 1 as ˆ H L g = 12 aL g ∆ ˆΦ (A9) ˆ H L J = ϕ I c a (cid:32) − cos (cid:32) ∆ ˆΦ ϕ (cid:33)(cid:33) (A10) ˆ H C g = C g a ˆ V C g = 12 aC g ˆ Q C g (A11) ˆ H C J = C J a (cid:32) ∂ ∆ ˆΦ ∂t (cid:33) (A12) Appendix B: Inductance of the unit cell It is useful to define an effective inductance that takesinto account the parallel effect of the geometric inductance L g and the Josephson capacitance C J . Keeping in mind2that the impedance of an inductor L for the mode n is Z L = jω n L while the impedance of a capacitor C for thesame mode is Z c = 1 /jω n C , we can write Z L eff ,n = 1 Z L g + 1 Z C J jω n L eff ,n = 1 jω n L g + jωC J L eff ,n = 1 L g − ω n C J = 1 − ω n L g C J L g (B1)It is found that the dispersion coefficient of the n -thmode ( Λ n ) can be defined by the relation L eff ,n = L g − ω n L g C J ≡ Λ n L g (B2)It is now possible to compute the total inductance of theelementary cell by calculating the parallel of the effectiveinductance L eff ,n and the Josephson inductance ˆ L J L n = 1ˆ L J + 1 L eff ,n = ˆ L J L eff ,n ˆ L J + L eff ,n (B3)Hence using equations (B2) and (6) the unit cell induc-tance can be written as ˆ L n = L J L g (cid:16) ∆ ˆΦ /ϕ (cid:17) L J (cid:16) ∆ ˆΦ /ϕ (cid:17) (cid:18) L g L J0 sin ( ∆ˆΦ /ϕ )( ∆ˆΦ /ϕ ) + 1 − L g C J ω n (cid:19) = L g L g L J0 sin ( ∆ˆΦ /ϕ )( ∆ˆΦ /ϕ ) + 1 − L g C J ω n = Λ n L g n L g L J0 sin ( ∆ˆΦ /ϕ )( ∆ˆΦ /ϕ ) (B4) Appendix C: Nonlinear time-dependent flux operator The nonlinear time-dependent component of the flux dif-ference operator can be found through the constitutiveequation of an inductor (8), where a mode decomposi-tion has been performed. It has to be noticed that theinductance can be mode-dependent. The current throughthe unit cell can be calculated exploiting the telegrapher’sequation (5), where the voltage drop to ground comes from(4) and the cell inductance of the n -th mode is found fromequation (B4). Hence, the AC current passing through theline results to be ˆ I ( z, t ) = (cid:88) n sgn ( n ) (cid:115) (cid:126) ω n L n N (cid:16) ˆ a n e i ( k n z − ω n t ) + H . c . (cid:17) (C1) the nonlinear time-dependent flux operator δ ˆΦ is then δ ˆΦ = (cid:88) n sgn ( n ) (cid:115) (cid:126) ω n L n N ˆ L n (cid:16) ˆ a n e i ( k n z − ω n t ) + H . c . (cid:17) = (cid:88) n sgn ( n ) (cid:114) (cid:126) ω n N (cid:112) Λ n L g ·· (cid:32) n L g L J sin ∆Φ DC + δ ˆΦ ϕ ∆Φ DC + δ ˆΦ ϕ (cid:33) − / ·· (cid:16) ˆ a n e i ( k n z − ω n t ) + H . c . (cid:17) == (cid:88) n n L g L J sin ∆Φ DC + δ ˆΦ ϕ ∆Φ DC + δ ˆΦ ϕ − / δ ˆΦ (0) n (C2)where we have identified δ ˆΦ (0) n ≡ sgn ( n ) (cid:114) (cid:126) ω n N (cid:112) Λ n L g (cid:16) ˆ a n e i ( k n z − ω n t ) + H . c . (cid:17) = c n (cid:16) ˆ a n e i ( k n z − ω n t ) + H . c . (cid:17) (C3)with c n = sgn ( n ) (cid:114) (cid:126) ω n N (cid:112) L g Λ n Equation (C2) is a recursive relation that involves δ ˆΦ andcan be solved imposing a solution to the lowest perturba-tive order, hence making the substitution δ ˆΦ (cid:55)→ δ ˆΦ (0) = (cid:80) n δ ˆΦ (0) n on the right hand side, so that we get δ ˆΦ = (cid:88) n n L g L J sin ∆Φ DC + δ ˆΦ (0) ϕ ∆Φ DC + δ ˆΦ (0) ϕ − / δ ˆΦ (0) n (C4)By invoking the Taylor expansion of the square root intoequation (C4) for δ ˆΦ (0) (cid:28) ϕ , one obtains δ ˆΦ = (cid:88) n (cid:34) q ,n + q ,n (cid:16) δ ˆΦ (0) (cid:17) + q ,n (cid:16) δ ˆΦ (0) (cid:17) ++ q ,n (cid:16) δ ˆΦ (0) (cid:17) + O (cid:16) δ ˆΦ (0) (cid:17) (cid:35) δ ˆΦ (0) n (C5)We stress that the terms q ,n , q ,n , q ,n and q ,n are co-efficients of a Taylor expansion and result to be functionsof the external bias conditions (i.e., of the constant fluxdifference ∆Φ DC ). It’s worth noting here how the low-est perturbative order approach adopted in equation (C4)takes into account interactions of modes at the first order,that means a single multimodes interaction. While thepower expansion truncation up to the third order in equa-tion (C5) limits our model to the interaction of a singlemode ( δ ˆΦ (0) n ) with up to three modes. For a quantitativecomparison of the power expansion approach (C5) respectto the bare nonlinearity (C4) see Figure 7.3 Appendix D: Coupling Coefficients Defined p = cos (∆Φ DC /ϕ ) and p = sin (∆Φ DC /ϕ ) : χ = N (cid:126) (cid:20) I c ϕ (cid:18) − cos (cid:18) ∆Φ DC ϕ (cid:19)(cid:19) + ∆Φ L g (cid:21) (D1) χ ( n )1 = ω n (cid:32) L g Λ n (cid:34)(cid:18) I c p + ∆Φ DC L g (cid:19) q ,n + (cid:18) I c p ϕ + 1 L g + C J ∆ ω n (cid:19) q ,n (cid:35)(cid:33) (D2) χ ( n,l,m )3 = (cid:114) (cid:126) L N (cid:112) ω n Λ n ω l Λ l ω m Λ m (cid:34) (cid:18) I c p + ∆Φ DC L g (cid:19) q ,n + (cid:18) I c p ϕ + 1 L g (cid:19) q ,n q ,l + − I c p ϕ q ,n q ,l q ,m + C J q ,n q ,l ∆ ω n ∆ ω m,l + q ,n q ,l ∆ ω l ∆ ω n,m ] (cid:35) (D3) χ ( n,l,m,s )4 = (cid:126) L N (cid:112) ω n Λ n ω l Λ l ω m Λ m ω s Λ s (cid:34) (cid:18) I c p + ∆Φ DC L g (cid:19) q ,n + 12 (cid:18) I c p ϕ + 1 L g (cid:19) (2 q ,n q ,l + q ,n q ,l ) + − I c p ϕ q ,n q ,l q ,m − I c p ϕ q ,n q ,l q ,m q ,s ++ C J q ,n q ,l (∆ ω m ∆ ω s,l + ∆ ω n ∆ ω m,l ) + q ,n q ,l ∆ ω l ∆ ω m,n + q ,n q ,l ∆ ω n ∆ ω m,l ] (cid:35) (D4) Appendix E: Time evolution of the ladder operators:coupled mode equations Startig from equation (11) and (12) it is possible to workout the dynamics of the system in 3WM and 4WM regime.In 3WM regime, hence using (11) to compute the Heisen-berg equation it is obtained dˆ a p d t = i (cid:126) (cid:104) ˆ H { p , s , i } , ˆ a p (cid:105) = − i (cid:104) χ p1 ˆ A p + χ { p , s , i } ˆ a s ˆ a i (cid:105) (E1) dˆ a s d t = i (cid:126) (cid:104) ˆ H { p , s , i } , ˆ a s (cid:105) = − i (cid:104) χ s1 ˆ a s + χ { p , s , i } ˆ A p ˆ a † i (cid:105) (E2) dˆ a i d t = i (cid:126) (cid:104) ˆ H { p , s , i } , ˆ a i (cid:105) = − i (cid:104) χ i1 ˆ a i + χ { p , s , i } ˆ A p ˆ a † s (cid:105) (E3)4While in 4WM regime, through (12) dˆ a p d t = i (cid:126) (cid:104) ˆ H { p , s , j } , ˆ a p (cid:105) = − i (cid:104)(cid:16) ξ p + ξ pp + 2 ξ pp ˆ A † p ˆ A p ++ ξ ps ˆ a † s ˆ A p + ξ pj ˆ a † j ˆ a j (cid:17) ˆ A p + 2 χ { p , p , s , i } ˆ A † p ˆ a s ˆ a j (cid:105) (E4) dˆ a s d t = i (cid:126) (cid:104) ˆ H { p , s , i } , ˆ a s (cid:105) = − i (cid:104)(cid:16) ξ s + ξ ss + 2 ξ ss ˆ a † s ˆ a s + ξ ps ˆ A † p ˆ A p ++ ξ sj ˆ a † j ˆ a j (cid:17) ˆ a s + χ { p , p , s , i } ˆ A p ˆ A p ˆ a † j (cid:105) (E5) dˆ a i d t = i (cid:126) (cid:104) ˆ H { p , s , i } , ˆ a i (cid:105) = − i (cid:104)(cid:16) ξ j + ξ jj + 2 ξ jj ˆ a † i ˆ a j + ξ pj ˆ A † p ˆ A p ++ ξ sj ˆ a † s ˆ a s (cid:17) ˆ a i + χ { p , p , s , j } ˆ A p ˆ A p ˆ a † s (cid:105) (E6)The systems composed by equations (E1), (E2), (E3)and by equations (E4), (E5), (E6) are known as QuantumLangevin Equations (also Coupled Mode Equations in theclassical regime), and their solutions determine the timeevolution of the modes interacting into the TWJPA. Theundepleted pump approximation describes a regime wherethe signal and idler modes can be considered small com-pared to the pump mode. This approximation allows tosolve analytically the system of coupled differential equa-tions by substituting the ladder operator of the pump modewith its classical counterpart, having defined (cid:115) (cid:126) ω p C g N ˆ a p (cid:55)→ A p (E7)as the classical voltage amplitude of Eq. (4).In the 4M case, equation (E7) can be substituted into (E4) (cid:115) C g N (cid:126) ω p d A p d t = − i (cid:34) ( ξ p + ξ pp ) (cid:115) C g N (cid:126) ω p A p ++ 2 (cid:18) C g N (cid:126) ω p (cid:19) | A p | ξ pp A p + ξ ps (cid:115) C g N (cid:126) ω p A p ˆ a † s ˆ a s ++ ξ pj (cid:115) C g N (cid:126) ω p ˆ a † j ˆ a j A p + 2 χ { p , p , s , j } (cid:115) C g N (cid:126) ω p A ∗ p ˆ a s ˆ a j (cid:35) d A p d t = − i (cid:34) ( ξ p + ξ pp ) A p + 2 C g N (cid:126) ω p | A p | ξ pp A p ++ ξ ps A p ˆ a † s ˆ a s + ξ pj ˆ a † j ˆ a j A p + 2 χ p,p,s,j A ∗ p ˆ a s ˆ a j (cid:35) (E8) Keeping the leading terms in (E8) we get d A p d t ≈ − i (cid:20) ( ξ p + ξ pp ) A p + 2 C g N (cid:126) ω p | A p | ξ pp A p (cid:21) = − i (cid:18) ( ξ p + ξ pp ) + 2 C g N (cid:126) ω p | A p | ξ pp (cid:19) A p = − i Ψ p A p (E9)and by solving this latter, one can derive A p ( t ) = | A p , | e − i Ψ p t (E10)with Ψ p = ξ p + ξ pp + 2 ξ pp C g N (cid:126) ω p | A p | (E11) | A p , | is the voltage amplitude at t = 0 , the time inwhich the mode enters in the non-linear medium. For sakeof simplicity we have assumed the initial phase of A p equalto zero. Similarly, the time evolution for the signal andidler annihilation operators can be written as dˆ a s d t = − i (cid:34)(cid:16) ξ s + ξ ss + 2 ξ ss ˆ a † s ˆ a s + ξ ps (cid:18) C g N (cid:126) ω p (cid:19) | A p | ++ ξ si ˆ a † j ˆ a j (cid:17) ˆ a s + χ { p , p , s , j } (cid:18) C g N (cid:126) ω p (cid:19) A ˆ a † j (cid:35) ≈ − i (cid:34) (cid:18) ξ s + ξ ps (cid:18) C g N (cid:126) ω p (cid:19) | A p | (cid:19) ˆ a s ++ χ { p , p , s , j } (cid:18) C g N (cid:126) ω p (cid:19) A ˆ a † j (cid:35) = − i (cid:34) Ψ s ˆ a s + χ { p , p , s , j } (cid:18) C g N (cid:126) ω p (cid:19) A ˆ a † j (cid:35) (E12)with Ψ s = ξ s + ξ ps (cid:18) C g N (cid:126) ω p (cid:19) | A p | (E13)In the co-rotating frame dˆ a s d t = − iχ { p , p , s , j } (cid:18) C g N (cid:126) ω p (cid:19) A (cid:0) ˆ a CR j (cid:1) † e i (Ψ s +Ψ j ) t = − iχ | A p , | (cid:0) ˆ a CR j (cid:1) † e − i (2Ψ p − Ψ s − Ψ j ) t = − iχ | A p , | (cid:0) ˆ a CR j (cid:1) † e − i Ψ t (E14)where equation (E10) has been exploited, having Ψ = 2Ψ p − Ψ s − Ψ j (E15)and introducing χ = χ { p , p , s , j } C g N (cid:126) ω p (E16)5The 3WM system can be solved through the same pro-cedure starting from equation (E1) dˆ a p d t = i (cid:126) (cid:104) ˆ H { p , s , i } , ˆ a p (cid:105) d A p d t = − i (cid:104) χ p1 ˆ A p + χ { p , s , i } ˆ a s ˆ a i (cid:105) ≈ − iχ p1 A p (E17)whose solution is A p ( t ) = | A p , | e − iχ p1 t (E18)equation (E2) becomes dˆ a s d t = − i (cid:34) χ s1 ˆ a s + χ { p , s , i } (cid:115) C g N (cid:126) ω p | A p , | ˆ a † s e − iχ p1 t (cid:35) (E19)in the co-rotating frame dˆ a s d t = − iχ { p , s , i } (cid:115) C g N (cid:126) ω p | A p , | ˆ a † i e − i ( χ p1 − χ s1 − χ i1 ) t = − iχ | A p , | ˆ a † i e − i Ψ t (E20) dˆ a i d t = − iχ { p , s , i } (cid:115) C g N (cid:126) ω p | A p , | ˆ a † s e − i ( χ p1 − χ s1 − χ i1 ) t = − iχ | A p , | ˆ a † s e − i Ψ t (E21)with Ψ = χ p1 − χ s1 − χ i1 (E22)and χ = (cid:115) C g N (cid:126) ω p χ { p , s , i } (E23) Appendix F: Squeezing The correlation of the signal and idler photons resultsin a so-called squeezed output field of a TWJPA. One candefine the thermal photon number as N ( ω ) = (cid:88) n (cid:16) (cid:104) ˆ a † ω ˆ a ω n (cid:105) − (cid:104) ˆ a † ω (cid:105) (cid:104) ˆ a ω n (cid:105) (cid:17) (F1)and the squeezing parameter as M ( ω ) = (cid:88) n (cid:16) (cid:104) ˆ a ω ˆ a ω n (cid:105) − (cid:104) ˆ a ω (cid:105) (cid:104) ˆ a ω n (cid:105) (cid:17) (F2) From the definitions (25) and (26) one can compute therelation between the squeezing spectrum ( S ), the thermalphoton number and the squeezing parameter S ( ω ) = (cid:88) n (cid:104) ∆ ˆ Y θ ( ω )∆ ˆ Y θ ( ω n ) (cid:105) = 1 + 2 N ( ω ) − | M ( ω ) | (F3)For a vacuum input state, the number of thermal photonscan be easily calculated through equations (16) N ( ω ) = | v ( ω, t ) | (F4)Again using (16), the squeezing parameter for a vacuuminput state can be written as M ( ω ) = (cid:88) n (cid:16) (cid:104) ˆ a ω ˆ a Ω n (cid:105) − (cid:104) ˆ a ω (cid:105) (cid:104) ˆ a Ω n (cid:105) (cid:17) = (cid:88) n (cid:104) ˆ a ω ˆ a Ω n (cid:105) = (cid:88) n (cid:104) vac | (cid:16) u ( ω, t )ˆ a ω, + iv ( ω, t )ˆ a † ω (cid:48) , (cid:17) ·· (cid:16) u (Ω n , t )ˆ a Ω n, + iv (Ω n , t )ˆ a † Ω (cid:48) n, (cid:17) | vac (cid:105) e − i Ψ t = iu ( ω, t ) e − i Ψ t (cid:88) n v (Ω n , t ) (cid:104) vac | ˆ a ω, (cid:16) ˆ a Ω n , (cid:17) † | vac (cid:105) = iu ( ω, t ) v ( ω, t ) e − i Ψ t = (cid:16) Ψ χ | A p , | g sinh ( gt ) − i χ | A p , | g sinh ( gt ) cosh ( gt ) (cid:17) e − i Ψ t = | u ( ω, t ) v ( ω, t ) | e − i (cid:16) arctan (cid:16) g Ψ coth ( gt ) (cid:17) +Ψ t (cid:17) == | M ( ω ) | e iθ (F5)Where we exploited v ( ω (cid:48) ) = v ( ω ) and identified thesqueezing angle as θ = − (cid:16) arctan (cid:16) g Ψ coth gt (cid:17) + Ψ t (cid:17) (F6)Hence one can easily find the relation between M ( ω ) and N ( ω ) as | M ( ω ) | = | u ( ω, t ) v ( ω, t ) | = | u ( ω, t ) | | v ( ω, t ) | = (cid:16) | v ( ω, t ) | + 1 (cid:17) | v ( ω, t ) | = N ( ω )[ N ( ω ) + 1] (F7)that is the maximum allowed by the Heisenberg un-certainty principle and implies that the amplification isquantum limited [27].6 [1] R. H. Hadfield, “Single-photon detectors for optical quan-tum information applications,” Nature Photonics , vol. 3,no. 12, pp. 696–705, 2009.[2] C. M. Natarajan, M. G. Tanner, and R. H. Hadfield, “Su-perconducting nanowire single-photon detectors: physicsand applications,” Superconductor Science and Technology ,vol. 25, p. 063001, apr 2012.[3] J. Ullom, W. Doriese, D. Fischer, J. Fowler, G. Hilton,C. Jaye, C. Reintsema, D. Swetz, and D. Schmidt,“Transition-edge sensor microcalorimeters for x-ray beam-line science,” Synchrotron Radiation News , vol. 27, no. 4,pp. 24–27, 2014.[4] D. Fukuda, G. Fujii, T. Numata, K. Amemiya,A. Yoshizawa, H. Tsuchida, H. Fujino, H. Ishii, T. Itatani,S. Inoue, and T. Zama, “Titanium-based transition-edgephoton number resolving detector with 98% detection effi-ciency with index-matched small-gap fiber coupling,” Op-tics Express , vol. 19, pp. 870–875, Jan 2011.[5] A. J. Miller, S. W. Nam, J. M. Martinis, and A. V.Sergienko, “Demonstration of a low-noise near-infraredphoton counter with multiphoton discrimination,” AppliedPhysics Letters , vol. 83, no. 4, pp. 791–793, 2003.[6] A. E. Lita, A. J. Miller, and S. W. Nam, “Counting near-infrared single-photons with 95% efficiency,” Optics Ex-press , vol. 16, pp. 3032–3040, Mar 2008.[7] R. Yang and H. Deng, “Fabrication of the impedance-matched Josephson parametric amplifier and the study ofthe gain profile,” IEEE Transactions on Applied Supercon-ductivity , pp. 1–1, 2020.[8] M. Malnou, D. A. Palken, L. R. Vale, G. C. Hilton, andK. W. Lehnert, “Optimal operation of a Josephson para-metric amplifier for vacuum squeezing,” Physical ReviewApplied , vol. 9, p. 044023, Apr 2018.[9] M. A. Castellanos-Beltran, K. D. Irwin, L. R. Vale, G. C.Hilton, and K. W. Lehnert, “Bandwidth and dynamicrange of a widely tunable Josephson parametric amplifier,” IEEE Transactions on Applied Superconductivity , vol. 19,pp. 944–947, June 2009.[10] C. Eichler and A. Wallraff, “Controlling the dynamic rangeof a Josephson parametric amplifier,” EPJ Quantum Tech-nology , vol. 1, no. 1, p. 2, 2014.[11] T. C. White, J. Y. Mutus, I.-C. Hoi, R. Barends, B. Camp-bell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, E. Jef-frey, J. Kelly, A. Megrant, C. Neill, P. J. J. O’Malley,P. Roushan, D. Sank, A. Vainsencher, J. Wenner, S. Chaud-huri, J. Gao, and J. M. Martinis, “Traveling wave paramet-ric amplifier with Josephson junctions using minimal res-onator phase matching,” Applied Physics Letters , vol. 106,no. 24, p. 242601, 2015.[12] C. Macklin, K. O’Brien, D. Hover, M. E. Schwartz,V. Bolkhovsky, X. Zhang, W. D. Oliver, and I. Siddiqi,“A near–quantum-limited Josephson traveling-wave para-metric amplifier,” Science , vol. 350, no. 6258, pp. 307–310,2015.[13] N. Zobrist, B. H. Eom, P. Day, B. A. Mazin, S. R.Meeker, B. Bumble, H. G. LeDuc, G. Coiffard, P. Szypryt,N. Fruitwala, I. Lipartito, and C. Bockstiegel, “Wide-bandparametric amplifier readout and resolution of optical mi- crowave kinetic inductance detectors,” Applied Physics Let-ters , vol. 115, no. 4, p. 042601, 2019.[14] S. Chaudhuri, D. Li, K. D. Irwin, C. Bockstiegel, J. Hub-mayr, J. N. Ullom, M. R. Vissers, and J. Gao, “Broad-band parametric amplifiers based on nonlinear kinetic in-ductance artificial transmission lines,” Applied Physics Let-ters , vol. 110, no. 15, p. 152601, 2017.[15] M. R. Vissers, R. P. Erickson, H.-S. Ku, L. Vale, X. Wu,G. C. Hilton, and D. P. Pappas, “Low-noise kinetic induc-tance traveling-wave amplifier using three-wave mixing,” Applied Physics Letters , vol. 108, no. 1, p. 012601, 2016.[16] A. B. Zorin, “Josephson traveling-wave parametric am-plifier with three-wave mixing,” Physical Review Applied ,vol. 6, p. 034006, Sep 2016.[17] A. B. Zorin, M. Khabipov, J. Dietel, and R. Dolata,“Traveling-wave parametric amplifier based on three-wavemixing in a Josephson metamaterial,” in ,pp. 1–3, June 2017.[18] T. M. Buehler, D. J. Reilly, R. P. Starrett, A. D. Greentree,A. R. Hamilton, A. S. Dzurak, and R. G. Clark, “Single-shot readout with the radio-frequency single-electron tran-sistor in the presence of charge noise,” Applied Physics Let-ters , vol. 86, no. 14, p. 143117, 2005.[19] A. Aassime, G. Johansson, G. Wendin, R. J. Schoelkopf,and P. Delsing, “Radio-frequency single-electron transistoras readout device for qubits: Charge sensitivity and backac-tion,” Physical Review Letter , vol. 86, pp. 3376–3379, Apr2001.[20] S. W. Henderson, Z. Ahmed, J. Austermann, D. Becker,D. A. Bennett, D. Brown, S. Chaudhuri, H.-M. S. Cho,J. M. D’Ewart, B. Dober, S. M. Duff, J. E. Dusatko,S. Fatigoni, J. C. Frisch, J. D. Gard, M. Halpern, G. C.Hilton, J. Hubmayr, K. D. Irwin, E. D. Karpel, S. S. Ker-nasovskiy, S. E. Kuenstner, C.-L. Kuo, D. Li, J. A. B.Mates, C. D. Reintsema, S. R. Smith, J. Ullom, L. R.Vale, D. D. V. Winkle, M. Vissers, and C. Yu, “Highly-multiplexed microwave SQUID readout using the SLACMicroresonator Radio Frequency (SMuRF) electronics forfuture CMB and sub-millimeter surveys,” in Millimeter,Submillimeter, and Far-Infrared Detectors and Instrumen-tation for Astronomy IX (J. Zmuidzinas and J.-R. Gao,eds.), vol. 10708, pp. 170 – 185, International Society forOptics and Photonics, SPIE, 2018.[21] X. Guo, C.-l. Zou, C. Schuck, H. Jung, R. Cheng, and H. X.Tang, “Parametric down-conversion photon-pair source ona nanophotonic chip,” Light: Science & Applications , vol. 6,no. 5, pp. e16249–e16249, 2017.[22] T. H. A. van der Reep, “Mesoscopic hamiltonian for joseph-son traveling-wave parametric amplifiers,” Physical ReviewA , vol. 99, p. 063838, Jun 2019.[23] A. B. Arne L. Grimsmo, “Squeezing and quantum stateengineering with Josephson travelling wave amplifiers,” npjQuantum Information , 2017.[24] R. Loudon, The Quantum Theory of Light, third edition .Oxford Science Publication, 2000.[25] U. Vool and M. Devoret, “Introduction to quantum electro-magnetic circuits,” International Journal of Circuit Theory and Applications , vol. 45, no. 7, pp. 897–934, 2017.[26] R. Gambini, “Parametric amplification with a trilinearHamiltonian,” Physical Review A , vol. 15, no. 3, pp. 1157–1168, 1977.[27] E. Jeffrey, D. Sank, J. Y. Mutus, T. C. White, J. Kelly,R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth,A. Megrant, P. J. J. O’Malley, C. Neill, P. Roushan,A. Vainsencher, J. Wenner, A. N. Cleland, and J. M. Marti-nis, “Fast accurate state measurement with superconduct-ing qubits,” Physical Review Letter , vol. 112, p. 190504,May 2014.[28] S. M. Barnett and P. M. Radmore, Methods in TheoreticalQuantum Optics . Clarendon Press - Oxford, 1997. ACKNOWLEDGMENTS This work has been partially funded by the SU-PERGALAX project in the framework of the H2020-FETOPEN-2018-2020 call and the Joint Research ProjectPARAWAVE of the European Metrology Programme forInnovation and Research (EMPIR). This project has re-ceived funding from the EMPIR programme co-financedby the Participating States and from the European UnionsHorizon 2020 research and innovation programme. AUTHOR CONTRIBUTION A.G., L.F., E.E. and L.C. conceptualized the work.A.G., L.F. and E.E. carried out the theoretical and numer-ical analysis. A.G., L.F. and E.E. wrote the manuscript.L.C. and A.M. participated in the discussion and editingof the manuscript.