Networks of nonlinear superconducting transmission line resonators

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 a 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, that 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 us to measure spatial correlations or perform manipulations for the 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–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–7, 10–13]. In this context, the nonlinearity [8, 9] can be regarded as an interaction between the 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. Since polaritons are joint excitations of two degrees of freedom, the cavity and the nonlinear quantum system, typically two species of polaritons emerge, which are separated in energy by the coupling between the cavity and the nonlinear quantum system. For the existence and addressability of the polaritons, it is therefore necessary to operate 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 analogue of cavity QED, superconducting circuit QED. In circuit QED microwave photons confined by quasi-one-dimensional (1D) transmission line resonators are coupled to Josephson-junction-based nonlinear quantum systems (qubits) [16–18]. The very small mode volume of the 1D transmission line and the very large dipole moment of the Josephson qubit result in exceptionally large coupling rates [27], whereas the superconducting gap ensures low dissipation.

In this work, we propose an even simpler approach, that 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]. The differences between the eigenfrequencies of the combined resonator–Josephson junction system are here comparable to the fundamental mode frequency, which enables a clear separation between the eigenmodes. As a consequence, neither the nonlinearity of the modes nor a photon tunneling between the adjacent resonators of a network lead to appreciable coupling between different eigenmodes. This leads to a system where the excitations of one eigenmode are dynamically well separated from those of other eigenmodes and can serve as a reliable, high precision quantum simulator.

Transmission lines that are intersected by Josephson junctions [28] are currently being used for measurement [19] or amplification [29, 30] purposes in the strong driving limit. Here, 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. With present-day circuit QED technology it has been demonstrated that hundreds of resonators can be fitted on a single chip [31] and that disorder of resonator frequencies can be kept at a tolerable level for tens of resonators [32]. Moreover, due to the size of the resonators, the Josephson junctions are spatially well separated and can be wired up individually, which enables tuning each junction independently. In addition to the typical many-body features, the resonator network is predestined for studying the effects of decoherence and driving [13, 36–39] 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 [33–35].

The paper is organized as follows. First, we present two versions of 1D 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.

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 1D 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 the adjacent resonators. In this way, 2D lattices can also 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 [14, 15]. The advantage of this coupling scheme is that each transmission line resonator can be probed individually. This advantage however only holds for the 1D networks.

Zoom In Zoom Out Reset image size

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],

$$\begin{eqnarray*} &&\sum\limits_{\left\langle i,j\right\rangle}J_{i,j}\left(a_i a_j^{\dagger}+a_i^{\dagger}a_j\right), \end{eqnarray*}$$

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 such as the respective mode, the coupling capacitance and the coupling inductance. Based on the 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,

$$\begin{equation} H_{\mathrm{BH}}=\sum\limits_i\left[ \omega a_i^{\dagger}a_i+ U a_i^{\dagger}a_i^{\dagger}a_ia_i+J\left(a_i^{\dagger}a_{i+1}+a_ia_{i+1}^{\dagger}\right)\right]. \end{equation} \tag{ 1 }$$

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.

3.1. 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 the 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 that 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,

$$\begin{eqnarray*} &&\mathcal{L} = \mathcal{L}_{\mathrm{l}} + \mathcal{L}_{\mathrm{r}} + \mathcal{L}_{\mathrm{JJ}} \end{eqnarray*}$$

where

$$\begin{eqnarray*} \mathcal{L}_{\mathrm{l/r}}&=&\sum_{i=1}^N\frac{\Delta \cdot c}{2}\dot{\phi}_{i,\mathrm{l/r}}^2-\sum_{i=1}^{N-1}\frac{(\phi_{i,\mathrm{l/r}}-\phi_{i+1,\mathrm{l/r}})^2}{2\Delta\cdot l}\\ \mathcal{L}_{\mathrm{JJ}}&=&\frac{C_{\mathrm{s}}}{2}\left(\dot{\phi}_{N,\mathrm{l}}-\dot{\phi}_{N,\mathrm{r}}\right)^2+E_{\mathrm{J}}\cos\left(\frac{2\pi}{\Phi_0}(\phi_{N,\mathrm{l}}-\phi_{N,\mathrm{r}})\right), \end{eqnarray*}$$

with Φ₀ = h/2e the quantum of flux, E_J the Josephson energy and $\phi _{i,\mathrm {l/r}}=\int _{-\infty }^t V_{i,\mathrm {l/r}}(t')\mathrm {d}t'$ the node fluxes [26]. $\mathcal {L}_{\mathrm {l/r}}$ is the Lagrangian of the transmission lines to the left and the right of the Josephson junction respectively. This is 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. $\mathcal {L}_{\mathrm {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.

Zoom In Zoom Out Reset image size

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 dimensions, 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 $\mathcal {L}_{\mathrm {l}}$ and $\mathcal {L}_{\mathrm {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,

$$\begin{eqnarray*} &&\mathcal{L}_{\mathrm{JJ}} = \mathcal{L}_{\mathrm{JJ}}^{\rm lin} + \mathcal{L}_{\mathrm{JJ}}^{\rm nonlin}, \end{eqnarray*}\noindent$$

where,

$$\begin{eqnarray*} \hphantom{,,}\mathcal{L}_{\mathrm{JJ}}^{\rm lin}&=&\frac{C_{\mathrm{s}}}{2}\left(\dot{\phi}_{N,{\rm l}}-\dot{\phi}_{N,{\rm r}}\right)^2-\frac{1}{2L_{\mathrm{J}}}\left(\phi_{N,{\rm l}}-\phi_{N,{\rm r}}\right)^2,\\ \mathcal{L}_{\mathrm{JJ}}^{\rm nonlin}&=&E_{\mathrm{J}}\cos\left(\frac{2\pi}{\Phi_0}(\phi_{N,{\rm l}}-\phi_{N,{\rm r}})\right)+\frac{1}{2L_{\mathrm{J}}}\left(\phi_{N,{\rm l}}-\phi_{N,{\rm r}}\right)^2\\ &=&E_{\mathrm{J}}\left(1+\sum\limits_{n=2}^{\infty}\frac{-1^n}{(2n)!}\left(\frac{2\pi(\phi_{N,{\rm l}}-\phi_{N,{\rm r}})}{\Phi_0}\right)^{2n}\right). \end{eqnarray*}\noindent$$

Here, L_J = Φ²₀/(4π²E_J) 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,

$$\begin{equation} \frac{-(\phi_{1,\mathrm{l/r}}-\phi_{2,\mathrm{l/r}})}{\Delta\cdot l} = \Delta \cdot c \ddot{ \phi}_{\mathrm{l/r},1}, \end{equation}\noindent \tag{ 2 }$$

the equations for the node fluxes at the Josephson junction are,

$$\begin{equation} \frac{(\phi_{N-1,\mathrm{l/r}}-\phi_{N,\mathrm{l/r}})}{\Delta\cdot l}\pm\frac{1}{L_{\mathrm{J}}}(\phi_{N,{\rm l}}-\phi_{N,{\rm r}})=\Delta\cdot c \ddot{\phi}_{N,\mathrm{l/r}}\pm C_{\mathrm{s}}(\ddot{\phi}_{N,{\rm l}}-\ddot{\phi}_{N,{\rm r}}), \end{equation}\noindent \tag{ 3 }$$

and for every remaining node flux we get,

$$\begin{equation} \frac{\phi_{i-1,\mathrm{l/r}}-2\phi_{i,\mathrm{l/r}}+\phi_{i+1,\mathrm{l/r}}}{\Delta^2}v^2=\ddot{\phi}_{i,\mathrm{l/r}}, \end{equation}\noindent \tag{ 4 }$$

where $v=1/\sqrt {lc}$ is the phase velocity in the transmission line.

3.2. 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, $\mathcal {L}=\sum (\cdots )\Delta \to \mathcal {L}=\int (\cdots )\mathrm {d}x$. Accordingly, differences between the 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,

$$\begin{equation} \left.-\partial_x\phi_{\mathrm{l/r}}(x)\right |_{x=\mp\,L/2}=0, \end{equation} \tag{ 5 }$$

$$\begin{equation} \left.-\partial_x\phi_{\mathrm{l}}(x)\right |_{x=x_{\mathrm{J}}}=-\left.\partial_x\phi_{\mathrm{r}}(x)\right |_{x=x_{\mathrm{J}}}=lC_{\mathrm{s}}\delta\ddot{\phi}+\frac{l}{L_{\mathrm{J}}}\delta\phi. \end{equation}\noindent \tag{ 6 }$$

The first boundary condition (5) is a consequence of charge conservation. The current, which is the spatial derivative of the flux, has to vanish at both ends of the transmission line. The second boundary condition (6) again states the fact of charge conservation, imposing that the current flowing into the Josephson junction at point x_J on the right-hand side has to exit the Josephson junction on the left-hand side. This current that flows through the Josephson junction is in turn related to the phase drop δϕ = (ϕ_l|_x=xJ − ϕ_r|_x=xJ) across the junction. The equations of motion for the node fluxes in the remaining nodes of the transmission line resonator yield the wave equation,

$$\begin{equation} \partial_x^2\phi_{\mathrm{l/r}}v^2=\partial_t^2\phi_{\mathrm{l/r}}, \end{equation} \tag{ 7 }$$

for the flux ϕ_l/r.

3.3. 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,

$$\begin{equation} \phi_{\mathrm{l/r}}(x,t)=f_{\mathrm{l/r}}(x)g(t), \end{equation} \tag{ 8 }$$

and derive the appropriate spatial normal modes that fulfill wave equation (7), in combination with the boundary conditions (5) and (6). Wave equation (7) leads to the following equations for f_l/r and g:

$$\begin{eqnarray*} &&f_{\mathrm{l/r}}''(x)+k^2f_{\mathrm{l/r}}(x)=0,\\ &&\ddot{g}(t)+\omega^2\,g(t) =0, \end{eqnarray*}$$

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,

$$\begin{equation} f_{\mathrm{l}}(x)=A_{\mathrm{l}}\cos\left(k\left(x+\frac{L}{2}\right)\right), \end{equation} \tag{ 9 }$$

$$\begin{equation} f_{\mathrm{r}}(x)=A_{\mathrm{r}}\cos\left(k\left(x-\frac{L}{2}\right)\right). \end{equation} \tag{ 10 }$$

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,

$$\begin{eqnarray*} &&\left.\partial_x f_{\mathrm{l}}(x)\right|_{x=0}=-A_{\mathrm{l}}k\sin\left(k\frac{L}{2}\right)=A_{\mathrm{r}}k\sin\left(k\frac{L}{2}\right)=\left.\partial f_{\mathrm{r}}(x)\right|_{x=0}. \end{eqnarray*}$$

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 the 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,

$$\begin{eqnarray*} &&-f_{\mathrm{l}}'(0)g(t)=-lC_{\mathrm{s}}\omega^2 (f_{\mathrm{l}}(0)-f_{\mathrm{r}}(0))g(t)+\frac{l}{L_{\mathrm{J}}}(f_{\mathrm{l}}(0)-f_{\mathrm{r}}(0))g(t). \end{eqnarray*}$$

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,

$$\begin{equation} \frac{\omega_n}{v}=\cot\left(\frac{\omega_n}{v}\frac{L}{2}\right)\,\frac{2l}{L_{\mathrm{J}}}\,\left(1-\frac{\omega_n^2}{\omega_p^2}\right), \end{equation} \tag{ 11 }$$

where $\omega _p = 1/\sqrt {L_{\mathrm {J}}C_{\mathrm {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,

$$\begin{equation} \int_{-L/2}^{L/2}f_n(x)f_m(x)=\frac{2\left(k_m\cos\left(k_n\frac{L}{2}\right)\sin\left(k_m\frac{L}{2}\right)-k_n\cos\left(k_m\frac{L}{2}\right)\sin\left(k_n\frac{L}{2}\right)\right)}{k_m^2-k_n^2}. \end{equation} \tag{ 12 }$$

Using the transcendental equation for the eigenfrequencies (11) we get,

$$\begin{eqnarray*} &&\frac{2\left(k_m\cos(k_n\frac{L}{2})\sin(k_m\frac{L}{2})-k_n\cos(k_m\frac{L}{2})\sin(k_n\frac{L}{2})\right)}{k_m^2-k_n^2}=-C_{\mathrm{s}}\delta f_n\delta f_m, \end{eqnarray*}$$

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₂ scalar product to account for the discontinuity at the Josephson junction,

$$\begin{equation} \begin{array}{l} c\displaystyle\int_{-L/2}^{L/2}f_n(x)f_m(x)\mathrm{d}x+C_{\mathrm{s}}\delta f_n \delta f_m = \delta_{n,m} \eta_n,\\ \eta_n=c\left(\displaystyle\frac{L}{2}+\displaystyle\frac{\delta f_n^2}{k_n^2}\displaystyle\frac{l}{2L_{\mathrm{J}}}\left(1+\displaystyle\frac{\omega_n^2}{\omega_p^2}\right)\right). \end{array} \end{equation} \tag{ 13 }$$

Here, we opted not to normalize the spatial eigenmodes but rather regard them as modes with different effective mass η_n. Now, we can reformulate the flux field in terms of the eigenmodes f_n,

$$\begin{eqnarray*} &&\phi(x,t)= \sum_{n=1}^{\infty}g_n(t)f_n(x), \end{eqnarray*}$$

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,

$$\begin{eqnarray*} &&\mathcal{L}=\sum_{n=1}^{\infty}\left(\frac{\eta_n}{2}\dot{g}_n^2-\frac{1}{2}\eta_n\omega_n^2 g_n^2\right)+\left(E_{\mathrm{J}}\cos\left(\frac{2\pi \delta\phi}{\Phi_0}\right)+\frac{1}{2L_{\mathrm{J}}}\delta \phi^2\right), \end{eqnarray*}$$

where $\delta \phi =\sum _n g_n(t)\delta 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,

$$\begin{eqnarray*} H&=&\sum_{n=1}^{\infty}\left(\frac{\pi_n^2}{2\eta_n}+\frac{1}{2}\eta_n \omega_n^2g_n^2\right)-\left(E_{\mathrm{J}}\cos\left(\frac{2\pi\delta\phi}{\Phi_0}\right)+\frac{1}{2L_{\mathrm{J}}}\delta\phi^2\right),\\ \pi_n&=&\eta_n\dot{\phi}_n. \end{eqnarray*}$$

We quantize the theory in the usual way by imposing canonical commutation relations $[\hat {\pi }_n,\hat {g}_m]=-i\delta _{n,m}$. The generalized coordinates $\hat {g}_n$ and momenta $\hat {\pi }_n$ can be expressed in terms of annihilation and creation operators a_n and a^†_n via the relations,

$$\begin{eqnarray*} \hat{\pi}_n&=&-\mathrm{i}\sqrt{\frac{\eta_n\omega_n}{2}}(a_n-a_n^{\dagger}),\\ \hat{g}_n&=&\frac{1}{\sqrt{2\eta_n\omega_n}}(a_n+a_n^{\dagger}), \end{eqnarray*}$$

with $[a_n,a_m^{\dagger }]=\delta _{n,m}$. Finally, we arrive at the following Hamilton operator:

$$\begin{equation} \begin{array}{l} H=\displaystyle\sum\limits_{n=1}^{\infty}\omega_n\left(a_n^{\dagger}a_n+\displaystyle\frac{1}{2}\right)+H^{\mathrm{nonlin}},\\ H^{\mathrm{nonlin}}=-E_{\mathrm{J}}(\cos(\delta\tilde{\phi})+\textstyle\frac{1}{2}\delta\tilde{\phi}^2), \end{array} \end{equation} \tag{ 14 }$$

with

$$\begin{eqnarray*} &&\delta\tilde{\phi}=\sum_{n=1}^{\infty}\lambda_n\left(a_n+a_n^{\dagger}\right)\quad \mathrm{and} \quad \lambda_n=\frac{2\pi\delta f_n}{\Phi_0\sqrt{2\eta_n \omega_n}}. \end{eqnarray*}\noindent$$

3.4. 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_ia^†_j + a^†_ia_j for i ≠ j, because we already diagonalized the complete linear part of the Hamiltonian. Secondly, density–density interactions of the type a^†_ia_ia^†_ja_j can be neglected because we assume all the modes except 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:

$$\begin{eqnarray*} &&\fl \cos\left(\sum\limits_{n=1}^{\infty} \lambda_n(a_n+a_n^{\dagger})\right)=\cos\left(\lambda_1(a_1+a_1^{\dagger})\right)\cos\left(\sum\limits_{n=2}^{\infty} \lambda_n(a_n+a_n^{\dagger})\right)\\ &&-\sin\left(\lambda_1(a_1+a_1^{\dagger})\right)\sin\left(\sum\limits_{n=2}^{\infty} \lambda_n(a_n+a_n^{\dagger})\right). \end{eqnarray*}$$

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,

$$\begin{eqnarray*} &&\cos\left(\sum_n \lambda_n (a_n+a_n^{\dagger})\right)\to \prod_n \cos\left(\lambda_n(a_n+a_n^{\dagger})\right). \end{eqnarray*}$$

The cosine terms can further be approximated by a sum of normally ordered operator products with equal numbers of annihilation and creation operators,

$$\begin{equation} \cos\left(\lambda_n (a_n+a_n^{\dagger})\right)\to \alpha_{n,0}+\alpha_{n,1} a_n^{\dagger} a_n+\alpha_{n,2} a_n^{\dagger}a_n^{\dagger} a_n a_n+\cdots. \end{equation} \tag{ 15 }$$

The expansion coefficients α_n,k in equation (15) are not the usual expansion coefficients of the cosine since one has to also take into account that 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],

$$\begin{eqnarray*} &&(a+a^{\dagger})^m=\sum\limits_{k=0}^{[m/2]}\sum\limits_{i=0}^{m-2k}\frac{m!(a^{\dagger})^i a^{m-2k-i}}{2^k k! i! (m-2k-i)!}, \end{eqnarray*}$$

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,

$$\begin{eqnarray*} &&\cos\left(\lambda_n(a_n+a_n^{\dagger})\right)\to \mathrm{e}^{-\lambda_n^2/2}\left(1-\lambda_n^2 a_n^{\dagger}a_n+\frac{\lambda_n^4}{4} a_n^{\dagger} a_n^{\dagger}a_na_n+\cdots\right). \end{eqnarray*}$$

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 λ²_n〈a^†_na_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,

$$\begin{eqnarray*} &&H^{\mathrm{nonlin}}\to-E_J\left(\left(1-\lambda_1^2 a_1^{\dagger}a_1+\frac{\lambda_1^4}{4} a_1^{\dagger} a_1^{\dagger}a_1a_1\right)\prod_{n=1}^{\infty}\mathrm{e}^{-\lambda_n^2/2}+\lambda_1^2 a_1^{\dagger}a_1\right). \end{eqnarray*}$$

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:

$$\begin{eqnarray*} &&H_1=\left(\omega_1-\delta \omega\right) a_1^{\dagger}a_1-U_{1} a_1^{\dagger}a_1^{\dagger}a_1a_1, \end{eqnarray*}$$

with

$$\begin{eqnarray*} &&\delta \omega_{1} = E_{\mathrm{J}} \lambda_1^2\left(1-\prod\limits_{n=1}^{\infty} \mathrm{e}^{-\lambda_n^2/2}\right),\\ &&U_{1} = E_{\mathrm{J}} \frac{\lambda_1^4}{4}\prod\limits_{n=1}^{\infty} \mathrm{e}^{-\lambda_n^2/2}, \end{eqnarray*}$$

where δω₁ is a small renormalization of the fundamental mode due to the nonlinearity.

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₁. If the self-inductance of the superconducting quantum interference device (SQUID) is negligible, it can be modelled as a single Josephson junction whose critical current I_c = 2I₀|cos(Φ/(2Φ₀))| (cf [24]) can be tuned via the external flux Φ threading the SQUID loop.

Further discussion will be twofold. We first consider the experimental limitations that 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-the-art 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₀ , 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 critical currents, particularly focusing on the experimentally well-accessible range of 8 × 10⁻⁷ A < I_c < 4 × 10⁻⁶ A. A summary of all relevant parameters is given in table 1.

In order to solve the transcendental equation for the spatial mode frequencies (11) we first reformulate it in terms of the phase φ associated with the wave vector of the spatial eigenmode, φ = k(L/2),

$$\begin{equation} \frac{\varphi}{1-\frac{\varphi^2}{\varphi_{\mathrm{p}}^2}}\frac{L_{\mathrm{J}}}{l L}=\cot(\varphi), \end{equation} \tag{ 16 }$$

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 = φ_n2v/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,

$$\begin{eqnarray*} \varphi_{\mathrm{n}} &=& (2n+1)\frac{\pi}{2} \quad n \in \mathbb{N},\\ \omega_{\mathrm{n}} &=& (2n +1)\frac{\pi v}{L}\quad n \in \mathbb{N}. \end{eqnarray*}$$

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)\frac {\pi v}{L}\approx \omega _{\mathrm {p}}$. The latter arise from their mutual interaction.

To calculate the nonlinearity parameter for the fundamental mode, U₁, we need to evaluate the formally infinite product, $\prod \nolimits _{n=1}^{\infty }\mathrm {e}^{-\lambda _n^2/2}$. Here, the superconducting gap provides a natural cutoff 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 $\prod _{n=1}^{n_{\mathrm {cutoff}}}\mathrm {e}^{-\lambda _n^2/2}$. Here, n_cutoff is the largest of all n that fulfills the equality,

$$\begin{eqnarray*} &&\hbar\frac{\pi v}{L}(2n-1)<\Delta(T)\quad\mathrm{and}\quad\Delta(T)=3.52\,k_{\mathrm{B}} T_{\mathrm{c}} \sqrt{1-\frac{T}{T_{\mathrm{c}}}} , \end{eqnarray*}$$

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 K 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 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₁, where we performed a rotating wave approximation and assumed all other modes to be in the vacuum state. Here we can also see that the higher modes could be employed as onsite nonlinear modes 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⁻⁷ A < I_c < 4 × 10⁻⁶ A. This results in a tuning range for the nonlinearity U₁ of $U_{\min }/(2\pi \hbar )=10\,{\rm KHz}$ to $U_{\max }/(2\pi \hbar )=4\,{\rm MHz}$ and a frequency drag for the fundamental mode of $\omega _{\min }/(2\pi )=5.2\,{\rm GHz}$ to $\omega _{\max }/(2\pi )=8.5\,{\rm 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.

Zoom In Zoom Out Reset image size

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 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⁻⁸ A < I_c < 4 × 10⁻⁴ 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₁a₁a₁a^†₂. The prefactor for this term in the Hamiltonian is of the order, E_Jλ³₁λ₂. This should be compared to the frequency difference between three times the fundamental mode frequency ω₁ and the second odd mode frequency ω₂. 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.

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.

We acknowledge valuable discussions with J Bourassa and K F Wulschner. This work was part of the collaborative research center SFB 631 and the Emmy Noether project HA 5593/1-1, both funded by the German Research Foundation (DFG). FD, AM and RG further acknowledge support from the German excellence initiative 'Nanosystems Initiative Munich' (NIM), the EU Marie Curie initial training network CCQED and the EU Specific Targeted Research Program PROMISCE.