Networks of nonlinear superconducting transmission line resonators

We investigate a network of coupled superconducting transmission line resonators, each of them made nonlinear with a capacitively shunted Josephson junction coupling to the odd flux modes of the resonator. The resulting eigenmode spectrum shows anticrossings between the plasma mode of the shunted junction and the odd resonator modes. Notably, we find that the combined device can inherit the complete nonlinearity of the junction, allowing for a description as a harmonic oscillator with a Kerr nonlinearity. Using a dc SQUID instead of a single junction, the nonlinearity can be tuned between 10 kHz and 4 MHz while maintaining resonance frequencies of a few gigahertz for realistic device parameters. An array of such nonlinear resonators can be considered a scalable superconducting quantum simulator for a Bose-Hubbard Hamiltonian. The device would be capable of accessing the strongly correlated regime and be particularly well suited for investigating quantum many-body dynamics of interacting particles under the influence of drive and dissipation.


Introduction
In present-day research the physics of quantum many-body systems has gained a lot of attention. Many-body Hamiltonians are investigated as minimal models for collective physical phenomena ranging from quantum phase transitions to quantum transport and nonequilibrium dynamics. Moreover, technological applications are being investigated for systems that can be described by interacting many-body Hamiltonians, such as topologically protected quantum states, novel schemes for quantum error correction or one way quantum computing.
There are only few situations where analytical solutions for interacting many-body Hamiltonians can be found. Normally one would use advanced numerical techniques to investigate theoretically the properties of the many-body Hamiltonian, but numerical methods are limited by the exponential growth of the Hilbert space of many-body systems as one increases the number of particles. As an alternative approach to circumvent this dilemma, quantum simulators are now being explored intensively. These are well controllable quantum systems that may emulate the physics of other systems which are less amenable to experimental investigation and are therefore expected to give rise to interesting and yet scarcely explored physics. Present day quantum simulators are therefore often strongly magnified versions of the respective system one is interested in, where for example lattice constants greatly exceed interatomic distances in a crystal. This can either allow to measure spatial correlations or perform manipulations for specific tasks such as quantum computing. Quantum many-body Hamiltonians can for example be simulated with cold atoms trapped by laser fields in various shapes and dimensions [1], in ion traps [2,3] or arrays of cavity quantum electrodynamics (QED) systems [4,5,6,7].
In the latter, instead of massive particles, joint excitations of the cavity field mode and a nonlinear quantum system, so-called polaritons, interact with each other [4,5,6,7,10,11,12,13]. In this context, the nonlinearity [9,8] can be regarded as an interaction between polaritons in the same cavity. Multiple cavities can be coupled to form a network, such that photons can hop between different cavities and consequently the polaritons are capable of propagating through the network of cavities. Polaritons however only form in the strong coupling regime, which can be reached by increasing the coupling between the nonlinear system and the cavity until it exceeds the relevant decay rates. This condition is nowadays extremely well met in the circuit analog of cavity QED, superconducting circuit QED. In circuit QED microwave photons confined by quasi one-dimensional transmission line resonators are coupled to Josephson junction based nonlinear quantum systems (qubits) [16,17,18]. The very small mode volume of the one-dimensional transmission line and the very large dipole moment of the Josephson qubit result in exceptionally large coupling rates [27] while the superconducting gap ensures low dissipation.
In this work, we propose an even simpler approach, which requires less control circuitry and therefore facilitates coupling multiple resonators to form scalable networks.
Instead of employing joint excitations of the resonator and the qubit, a transmission line resonator can directly be rendered nonlinear by intersecting it with a Josephson junction [19]. We show that the modes of such a device can achieve a significant amount of nonlinearity on the single photon level for realistic circuit parameters. If several of such transmission line resonators are connected to form a network, the dynamics of the photons propagating through this network is thus described by a Bose-Hubbard Hamiltonian. In this way our approach leads to quantum simulations of the physics of massive, interacting particles by employing massless photons. Furthermore, the setup would reduce the required amount of control circuitry and hence significantly increase the scalability of the setup. In addition to the typical many body features, the resonator network is predestined for studying the effects of decoherence and driving [13,31,32,33,34] as well as the transition from a single entity to a many body physics. In this sense, the system we propose can also be viewed as a quantum simulator. Notably, a set of experimental analysis methods for the quantum statistics of resonator networks have been demonstrated in recent years [28,29,30] The article is organized as follows. First we present two versions of one-dimensional arrays of coupled nonlinear transmission line resonators that can be described by a Bose-Hubbard Hamiltonian. In the next step we derive the full Hamiltonian for a single transmission line resonator intersected by a Josephson junction. We then focus our attention on the fundamental mode of the resonator, derive an approximate Hamiltonian for it that includes the junction nonlinearity and show that it represents a single site of a Bose-Hubbard Hamiltonian. Finally, regarding an experiment we quantify how much nonlinearity can be achieved with parameters of recently realized circuit QED setups.

Chains of transmission line resonators
Circuit QED offers ample possibilities to couple several resonators to a network as there are basically no limitations on the topology and geometry of networks of transmission line resonators [20,21], except for constraints imposed by the detection circuitry and by space on the chip. In figure 1 we display two specific configurations of one-dimensional networks (chains) of transmission line resonators. In figure 1 a), the resonators are coupled at both ends [25], which form a small capacitance and thereby enable photon hopping between adjacent resonators. In this way, also two-dimensional lattices can be formed by coupling more than two resonators at their ends [21]. The advantage of this coupling scheme is scalability whereas a drawback results from the fact that one can only probe the output fields at the borders of the entire network [25]. However, transmission line resonators can also be coupled by reducing the distance between their central conductors for a certain length as shown in figure 1 b). Depending on the exact position of the convergence, and the specific mode under consideration, the coupling between the two resonators is in general both capacitive and inductive [15,14]. The advantage of this coupling scheme is that each transmission line resonator can be probed individually. This advantage however only holds for one-dimensional networks. For both designs the coupling constants are typically smaller than the frequencies of the field modes and therefore a rotating wave approximation is applicable. The coupling Hamiltonian consequently reads [20], where the a i are the mode operators of the field modes in the respective transmission line resonator, the sum comprises all coupled transmission line resonators and J i,j incorporates all microscopic details of the interaction like the respective mode, the coupling capacitance and the coupling inductance. Based on present day precision in sample production we assume, for the following, all coupling constants to be equal, J i,j = J. As shown below every nonlinear transmission line resonator can be modeled by a harmonic oscillator Hamiltonian with a Kerr nonlinearity U. For the whole chain of nonlinear transmission line resonators we therefore get a Bose-Hubbard Hamiltonian of the form,

Nonlinear transmission line resonator
In this section, we derive the mode spectrum of a single transmission line resonator that is intersected by a Josephson junction [22]. To this end, we first disregard its nonlinearity and compute the normal modes of the harmonic part of its Hamiltonian. The nonlinearity is then written in terms of these normal modes. To correctly identify the boundary conditions and the wave equation governing the flux function of the nonlinear resonator, we start by deriving the discretized, lumped-element representation of the transmission line resonator.

Lumped-element model
We divide the resonator into finite length lumps that are represented in our circuit diagram by nodes. These nodes are inductively connected in a series. The inductance of each connection is ∆ · l, the length of the unit cell ∆ times the inductance per unit length l of the transmission line. Each lump of transmission line is itself capacitively coupled to the ground plane. We represent this by a capacitance to ground at each node in the circuit diagram. The capacitance ∆ · c is again the product of the length of the unit cell ∆ and the capacitance per unit length c. The Josephson junction interrupts the transmission line at the location x J , where the two parts of the transmission line are connected capacitively and through the overlap of their Cooper-pair wavefunctions. The latter is represented in the circuit diagram by a nonlinear Josephson inductance which is shunted by the capacitance of the Josephson junction itself and an additional shunting capacitance C s = C J + C shun . The Lagrangian of the whole discretized transmission line setup reads, where, [26]. L l/r are the Lagrangians of the transmission lines to the left and the right of the Josephson junction respectively. These are represented by N nodes each, with the first node being located at the respective end of the transmission line and the Nth node being the node next to the Josephson junction, see figure 2. L JJ is the Lagrangian of the Josephson junction. Our goal is to find the eigenmodes for the linear part of the Hamiltonian which will be quantized in the usual way.
We will proceed to accomplish this task as follows: We first separate the linear parts of the Langrangian from the rest, take the continuum limit and derive the wave equation for the linear part. Then we do a separation ansatz for the spatial and temporal dimension, identify the spatial boundary conditions and derive the corresponding eigenmodes. The time variation of the eigenmodes will then be described by harmonic oscillator equations of motion and therefore can be quantized in the same way as independent harmonic oscillators.
Since L l and L r are linear Langrangians, the separation of the linear and nonlinear parts of the Lagrangian can be accomplished by dividing the Lagrangian for the Josephson junction in the following way, is the Josephson inductance. The linear part of the Lagrangian leads to the following set of equations of motion for the node fluxes. For the nodes at both ends of the transmission line the equations of motion read, the equations for the node fluxes at the Josephson junction are, and for every remaining node flux we get, where v = 1/ √ lc is the phase velocity in the transmission line.

Continuum limit
To take the continuum limit, we shrink the size, ∆, of our unit-cell of transmission line and transform every sum in the Lagrangian into an integral over the spatial coordinate, L = (...)∆ → L = (...)dx. Accordingly, differences between neighboring node fluxes divided by the size of our unit cell become derivatives with respect to the spatial dimension. We note that the node fluxes to the right and the left of the Josephson junction can have a finite difference in the continuum limit which is a consequence of our model assumption that the Josephson junction has no spatial extension in the direction of the transmission line. This is a valid assumption as long as we study spatial modes of a wavelength that is much larger than the size of the Josephson junction. In this way of downsizing the unit cell, the inductance of the Josephson junction remains finite, even though its extension becomes vanishingly small. The equations of motion for the fluxes at discontinuities of the transmission line (3) and (2) provide us with the boundary conditions, The first boundary condition (5) for the flux φ l/r .

Eigenmodes of the linearized transmission line resonator
To find the eigenmodes of the linearized transmission line resonator we use an ansatz with separation of variables, and derive appropriate spatial normal modes that fulfill the wave equation (7), in combination with the boundary conditions (5) and (6). The wave equation (7) leads to the following equations for f l/r and g, where ω and k are the frequency and wave vector of the normal mode which are connected by the linear dispersion relation |k|v = ω. The boundary conditions at both ends of the transmission line resonator directly translate into boundary conditions for the function f (x). To ensure current conservation at both ends of the transmission line resonator (5) we choose the functions f l and f r to read, We now restrict our calculations to the case where the Josephson junction is in the middle of the transmission line resonator, x J = 0, and hence current conservation through the Josephson junction, compare equation (6), requires, There are two different ways to fulfill this boundary condition. Choosing k n = 2πn/L provides us with symmetric (even) spatial modes that do not couple to the Josephson junction because no current is flowing through it. Alternatively we can choose A l = −A r which provides us with the antisymmetric (odd) spatial modes. We specifically choose A l = −A r = 1 and consequently obtain orthogonal but not normalized modes. The remaining work to be done is to connect the current flowing through the Josephson junction and the phase drop across the Josephson junction with the help of boundary condition (6), where we insert our separation ansatz (8) and get, Here, we used the equation of motion for the function g(t). With our ansatz for the spatial modes, equations (9) and (10), we get the transcendental equation for the spatial mode frequencies ω n , where ω p = 1/ √ L J C s is the plasma frequency of the capacitively shunted Josephson junction. With the help of this transcendental equation we compute the frequencies for the odd modes of the system with n = 2k + 1. For the even modes the current through the Josephson junction vanishes and consequently the flux drop across the Josephson junction also vanishes.
One can check the orthogonality of the odd modes by integrating the product of two different (m = n) odd modes over the length of the transmission line resonator, Using the transcendental equation for the eigenfrequencies (11) we get, where δf n is the flux drop of the spatial eigenmode n, δf n = f n,l (0) − f n,r (0). This result leads us to alter the usual l 2 scalar product to account for the discontinuity at the Josephson junction, Here, we opted not to normalize the spatial eigenmodes but rather regard them as modes with different effective masses η n . Now we can reformulate the flux field in terms of the eigenmodes f n , and with the help of the orthogonality relation (13) we rearrange the integral in the Lagrange function into a sum over the infinite number of spatial eigenmodes, where δφ = n g n (t)δf n . Here, we already neglected the even spatial modes, because they are completely decoupled from the odd modes. The corresponding Hamilton function can be found by a Legendre transformation which leaves the nonlinear part of the Lagrangian unaltered, We quantize the theory in the usual way by imposing canonical commutation relations [π n ,ĝ m ] = −iδ n,m . The generalized coordinatesĝ n and momentaπ n can be expressed in terms of annihilation and creation operators a n and a † n via the relations, with a n , a † m = δ n,m . Finally, we arrive at the following Hamilton operator, ω n (a † n a n + with δφ = ∞ n=1 λ n a n + a † n and, λ n = 2πδf n Φ 0 √ 2η n ω n .

Derivation of the Kerr-nonlinearity
In principle, every odd mode is coupled to all other odd modes by the nonlinear Hamiltonian (14). As we illustrate in the following, the interactions between different modes are however negligible in the low energy limit. Firstly, we note that the nonlinear part of the Hamiltonian (14) does not contain any hopping terms of the form a i a † j + a † i a j for i = j because we already diagonalized the complete linear part of the Hamiltonian. Secondly, density-density interactions of the type a † i a i a † j a j can be neglected because we assume all the modes but the fundamental mode to be in the vacuum state and the negligible interaction between the modes prevents them from becoming populated. Finally, the remaining higher order terms do not contribute a significant coupling. To separate the Hamilton function of the fundamental mode from the rest, we apply a rotating wave approximation, the validity of which we confirm below. This approximation allows us to substantially simplify the cosine term of the nonlinear part of the Hamilton operator (14). To this end we, consider the identity, λ n (a n + a † n ) .
The sine terms will only contain operator products with odd powers of the individual mode operators, hence the second addend can be neglected as part of the rotating wave approximation. If we iterate this consideration, we can approximate the cosine of the sum of our mode position operators by the product of the cosines of the mode position operators, cos n λ n (a n + a † n ) → n cos λ n (a n + a † n ) .
The cosine terms can further be approximated by a sum of normally ordered operator products with equal numbers of annihilation and creation operators, cos λ n (a n + a † n ) → α n,0 + α n,1 a † n a n + α n,2 a † n a † n a n a n + . . . .
The expansion coefficients α n,k in equation (15) are not the usual expansion coefficients of the cosine since one has to take into account that also higher powers of λ n (a n + a † n ) generate operator products of lower powers if we attempt to normal order the mode operator binomials. With the formula for mode operator binomials given in [23], where [m/2] denotes the largest integer less than or equal to m/2, we can identify the prefactors for the rotating wave expansion of the cosine, cos λ n (a n + a † n ) → e − λ 2 n 2 1 − λ 2 n a † n a n + λ 4 n 4 a † n a † n a n a n + . . . .
We assume that only the fundamental mode is populated and can therefore neglect all modes with n = 1 in H nonlin . Since |λ n | ≪ 1, as will be confirmed below, and because we are only interested in low photon numbers where λ 2 n a † n a n ≪ 1, we can truncate the cosine expansion after the quartic order term. We thus get for the nonlinear part of our Hamilton operator in the rotating wave and low photon number approximation, Any constant in the Hamilton operator provides us with an unmeasurable shift of the overall energy, therefore we may omit all constant terms and get the following Hamiltonian for the fundamental mode in rotating wave approximation, where δω 1 is a small renormalization of the fundamental mode due to the nonlinearity.

Experimental prerequisites
In this section we discuss how one would realize the above introduced nonlinear transmission line resonator. We propose to use a dc SQUID instead of the single Josephson junction in order to have a tunable and large nonlinearity U 1 . If the self inductance of the SQUID is negligible, it can be modeled as a single Josephson junction whose critical current I c = 2I 0 | cos(Φ/(2Φ 0 ))| (cf. [24]) can be tuned via the external flux Φ threading the SQUID loop.
The further discussion will be twofold. We first consider the experimental limitations which are imposed on our theoretical model by the physics of the superconducting circuitry and then discuss how the setup should be tuned to realize a sufficient amount of nonlinearity. In an experiment, it is favorable that the fundamental mode of the nonlinear resonator is within a frequency band suitable for state-of-theart microwave measuring devices. Furthermore, the thermal population should be insignificant at least at the millikelvin temperatures used in most experiments. We therefore choose the length L of the transmission line resonator such that for all relevant critical currents I c the frequency of the fundamental mode is in the range 4 GHz < ω/(2π) < 8 GHz. All other design parameters of the transmission line (Z 0 , l, c) are typical values found in many circuit QED setups. The critical current I c of the SQUID can be varied in two different ways. Either by design, or, as sketched above, by applying a flux bias to the SQUID loop. In this work we study a wide range of In order to solve the transcendental equation for the spatial mode frequencies (11) we first reformulate it in terms of the phase ϕ associated to the wave vector of the spatial eigenmode, ϕ = k(L/2), where ϕ p is the corresponding phase for the capacitively shunted dc-SQUID, ϕ p = (ω p /v)(L/2). Every time the left hand side of this equation equals the right hand side for some specific ϕ = ϕ n , we have found one valid mode frequency ω n = ϕ n 2v/L. For our choice of parameters, the right hand side has a narrow divergence for ϕ = ϕ p respectively ω = ω p and quickly approaches zero for values of ϕ unequal to ϕ p . Consequently, we have two different possibilities for approximate solutions of the transcendental equation. For ϕ = ϕ p , we can approximate the left hand side of the transcendental equation by 0 and find, This corresponds to the situation where the plasma frequency of the capacitively shunted dc-SQUID ω p is off-resonant with the unperturbed frequency of the transmission line resonator mode and consequently we recover the unperturbed frequencies of the resonator without dc-SQUID. The situation changes if one searches for system modes near the plasma frequency of the capacitively shunted dc-SQUID, ϕ ≈ ϕ p . Then, the divergence of the left hand side becomes infinitely narrow and the solution of the transcendental equation is approximately ϕ n ≈ ϕ p . These considerations are confirmed by the numerically calculated mode spectrum of the linear part of the combined device.
Here, we find the spectrum of the odd modes of the unperturbed transmission line resonator modes and the plasma frequency of the capacitively shunted dc-SQUID with avoided crossings near (2n + 1) πv L ≈ ω p . The latter arise from their mutual interaction.  To calculate the nonlinearity parameter for the fundamental mode, U 1 , we need to evaluate the formally infinite product, ∞ n=1 e − λ 2 n 2 . Here the superconducting gap provides a natural cut-off frequency. If the first excited state of a mode of the transmission line resonator exceeds the superconducting gap in energy, Cooper pairs will break and strong dissipation processes will start. Therefore we only include modes that fit energetically into the superconducting gap in the above product and get n cutoff n=1 e − λ 2 n 2 . Here, n cutoff is the largest of all n that fulfills the equality, with ∆ the energy gap of the superconductor. Here we used the critical temperature of niobium. Niobium technique is favorable because of its high transition temperature of 9.2 Kelvin that enables fast precharacterization at liquid helium temperatures before cooling the whole device down to the millikelvin range. The nonlinearity of the fundamental mode is inherited from the capacitively shunted dc-SQUID for values of the critical current I c where the frequency of the fundamental mode is near the plasma frequency. This principle of inheritance, where a collective mode of the composite resonator-junction system adopts the properties of one of its constituents if its frequency matches the frequency of the respective isolated constituent, is further illustrated by the nonlinearity U n of higher modes (see figure 3). The nonlinearities for these modes are calculated in the same manner as the nonlinearity for the fundamental mode, U 1 , where we performed a rotating wave approximation and assumed all other modes to be in the vacuum state. Here we can see that also higher modes could be employed as onsite nonlinear mode as they acquire the same amount of nonlinearity. In figure 3, we hatched a proposed tuning range for the SQUID critical current of 8 · 10 −7 A < I c < 4 · 10 −6 A. This results in a tuning range for the nonlinearity U 1 of U min /(2π ) = 10 kHz to U max /(2π ) = 4 MHz and a frequency drag for the fundamental mode of ω min /(2π) = 5.2 GHz to ω max /(2π) = 8.5 GHz. This provides us with sufficient nonlinearity to exceed present day decay rates of transmission line resonators γ/(2π) ≈ 100 kHz while still keeping the frequency drag within the bandwidth of state of the art microwave detection.
These values allow one to investigate the Bose-Hubbard Hamiltonian of (1) with in situ tunability between the onsite nonlinearity dominated regime U > J and the hopping dominated regime J > U using present day technology. Finally, we note that the nonlinearity U max is here bound from above by the charging energy (see figure 3). Therefore, it could be further increased by decreasing the shunting capacitance C s . However, this would shift the anticrossing marking the region of tunability, and hence the operating interval, to lower critical currents, finally resulting in fabrication and noise issues.
In the whole presented range of critical currents 8 · 10 −8 A < I c < 4 · 10 −4 A we checked the validity of our rotating wave approximation. The most dominant process of all neglected ones is the exchange of three fundamental mode excitations with one excitation of the next higher odd mode, a 1 a 1 a 1 a † 2 . The prefactor for this term in the Hamiltonian is of the order, E J λ 3 1 λ 2 . This should be compared to the frequency difference between three times the fundamental mode frequency ω 1 and the second odd mode frequency ω 2 . For the whole considered range of critical currents the ratio of the prefactor and the frequency difference is vanishingly small, which validates our assumption of decoupled field modes.

Summary
In summary, we have introduced a scheme to study interacting quantum many-body systems formed by photons. Our approach employs a network of transmission line resonators which are themselves nonlinear. Hence, no qubit circuits are required to generate a nonlinearity for microwave photons. Our approach thus reduces the amount of control circuitry required for building the cavity network. This reduced complexity is expected to enhance the perspectives for scaling up the technology to large networks of coupled circuit cavities. In combination with significant nonlinearities tunable between 10 kHz and 4 MHz, our system is a promising candidate for building quantum simulators for interesting few-and many-body physics problems.