Ultrastrong coupling dynamics with a transmon qubit

Cavity quantum electrodynamics allows for the study of light–matter interaction at the level of single atoms interacting with a single photon, both confined in a high-quality cavity. In practice, this interaction is typically due to the coupling of the light's electric field to the electric dipole moment of the atom [1]. When only a single mode of light and only two atomic levels are relevant, this situation can be described by the Jaynes–Cummings Hamiltonian (${\hslash }=1$),

$$\begin{eqnarray}&&{H}_{\mathrm{JC}}={\omega }_{r}\,{a}^{\dagger }a+\displaystyle \frac{{\omega }_{a}}{2}{\sigma }_{z}+g({a}^{\dagger }{\sigma }_{-}+a{\sigma }_{+}),\end{eqnarray} \tag{ 1 }$$

where ${\omega }_{r}$ is the cavity frequency, ${\omega }_{a}$ the atomic frequency and g the electric-dipole coupling. In this expression, a (${a}^{\dagger }$) is the photon annihilation (creation) operator and ${\sigma }_{i}$ are the Pauli matrices for the atomic levels. The Jaynes–Cummings Hamiltonian describes the exchange of a single quanta between the field and the atom leading to Rabi oscillations with angular frequency $2g$. The strong coupling regime is achieved when the coupling, g, is much larger than the dissipation rates of the system. This Jaynes–Cummings Hamiltonian can be realized with a wide variety of physical systems such as Rydberg atoms [1], quantum dots [2], trapped ions [3, 4] and superconducting circuits [5].

The Jaynes–Cummings Hamiltonian is, however, only an approximation of the Rabi Hamiltonian describing coupling between the cavity electric field, ${E}_{0}({a}^{\dagger }+a)$, and the atomic dipole moment, ${d}_{0}{\sigma }_{x}$,

$$\begin{eqnarray}&&{H}_{\mathrm{Rabi}}={\omega }_{r}\,{a}^{\dagger }a+\displaystyle \frac{{\omega }_{a}}{2}{\sigma }_{z}+g({a}^{\dagger }+a){\sigma }_{x},\end{eqnarray} \tag{ 2 }$$

where $g={d}_{0}{E}_{0}$. The Jaynes–Cummings Hamiltonian is a good approximation to ${H}_{\mathrm{Rabi}}$ when the coupling, g, is smaller than the system frequencies, $g\ll {\omega }_{a},{\omega }_{r}$. In this situation, the fast rotating term ${a}^{\dagger }{\sigma }_{+}+a{\sigma }_{-}$ appearing in the Rabi Hamiltonian can safely be dropped using the rotating wave approximation (RWA) and we recover equation (1). While more challenging to realize, there has recently been much attention to the situation where this approximation is no longer valid. This so-called ultrastrong regime, realized when the coupling strength approaches the system frequencies, differs remarkably from the Jaynes–Cummings regime [6–13]. Most significantly, while the ground state of ${H}_{\mathrm{JC}}$ is simply the product of the atomic ground state and vacuum of the field, the ground state of the Rabi Hamiltonian is an entangled atom-field state with a non-zero average photon number. From a practical point of view, this regime could also be useful in the context of quantum information processing [14–17].

Superconducting quantum circuits form a promising platform to realize and study this novel light–matter coupling regime. In particular, the realization of the ultrastrong coupling regime with flux qubits, acting as the atom, coupled to a microwave resonator have been theoretically studied [18] and experimentally implemented [19–22]. These experiments have primarily probed the spectral properties of the ultrastrong regime. A next step is to probe the dynamics of the system in this regime and, moreover, to probe its non-trivial ground state. With the system starting in its ground state, an approach is to non-adiabatically tune the coupling strength g from the ultrastrong coupling regime to the strong coupling regime. The system will readjust to this change by emitting photons as it relaxes back to its new ground state. Observing these photons would constitute a clear signature of the non-trivial nature of the ultrastrong coupling ground state. With system frequencies around 10 GHz [18–20, 22], this however requires changes in system parameters of the order of 10 pico-seconds. In practice, this therefore appears to be extremely challenging.

In this work, we address this problem by working with a low-frequency mode of a microwave cavity. We focus on the transmon qubit [23] capacitively coupled to an array of Josephson junctions realizing an inductance with a dissipationless impedance larger than the resistance quantum [24–26]. Together with its capacitance to ground, this superinductance plays the role of a multi-mode cavity. With $g/{\omega }_{r}\propto \sqrt{{Z}_{0}}$, where Z₀ is the cavity impedance [18, 27], this approach allows for large qubit-mode coupling strengths. Moreover, by using a low-frequency mode of the array, it is possible to realize ultrastrong coupling with only a moderately large coupling strength. In this situation, fast changes of system parameters are possible and allow for the observation of signatures of the ultrastrong coupling in the dynamics of the combined system. These dynamics can then be probed by taking advantage of the presence of multiple modes of the array and their cross-Kerr interaction [24].

The paper is organized as follows: we derive in section 2 the Hamiltonian of the system, taking into account the multi-mode structure of the array. In section 3, we identify parameters to reach the ultrastrong coupling regime. In section 4 the dynamics of the ultrastrong coupling regime are investigated. Finally, section 5 concludes the paper.

We consider the circuit of figure 1 which consists of a transmon qubit [23] coupled to an array of N Josephson junctions [24]. The transmon qubit is characterized by the Josephson energy E_J and the capacitance C_s, which for simplicity we take to include both the shunt capacitance and the junction capacitance. We assume the junctions forming the array to have a large Josephson energy such that, to a good approximation, they behave as weakly nonlinear inductances. These junctions are then characterized by their Josephson inductance L_J and junction capacitance C_J. Following [28, 29], the nonlinearity of the array junctions will be perturbatively reintroduced at a later step. Moreover, we take into account the capacitance to ground C₀ of the islands formed between the array junctions. The capacitance C_q couples the qubit to the array and will largely control their interaction strength. Finally, C_i is a capacitance to an external control field which will be used to probe the system and C_e is the capacitance to ground of the last array island, which can be constructed arbitrarily [24].

Zoom In Zoom Out Reset image size

Following the standard approach [30], the circuit Lagrangian reads

$$\begin{eqnarray}{ \mathcal L } & = & \displaystyle \sum _{n=0}^{N-1}\left[\displaystyle \frac{{C}_{J}}{2}{({\dot{\phi }}_{n}-{\dot{\phi }}_{n+1})}^{2}-\displaystyle \frac{1}{2{L}_{J}}{({\phi }_{n}-{\phi }_{n+1})}^{2}+\displaystyle \frac{{C}_{0}}{2}{\dot{\phi }}_{n}^{2}\right]\\ & & +\displaystyle \frac{{C}_{i}}{2}{\dot{\phi }}_{0}^{2}+\displaystyle \frac{{C}_{e}}{2}{\dot{\phi }}_{N}^{2}+\displaystyle \frac{{C}_{q}}{2}{({\dot{\phi }}_{0}-{\dot{\phi }}_{{qb}})}^{2}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} & & +\displaystyle \frac{{C}_{s}}{2}{\dot{\phi }}_{{qb}}^{2}+{E}_{J}\cos ({\phi }_{{qb}}\,/{\varphi }_{0})\\ & = & {\dot{\vec{\phi }}}^{T}\displaystyle \frac{{\boldsymbol{C}}}{2}\dot{\vec{\phi }}-{\vec{\phi }}^{T}\displaystyle \frac{{\boldsymbol{L}}}{2}\vec{\phi }+{E}_{J}\cos ({\phi }_{{qb}}\,/{\varphi }_{0}),\end{eqnarray} \tag{ 4 }$$

where ${\phi }_{n}$ is the node flux of the nth island of the array, ${\phi }_{{qb}}$ the node flux of the transmon's island and ${\varphi }_{0}={{\rm{\Phi }}}_{0}/2\pi$ with the magnetic flux quantum ${{\rm{\Phi }}}_{0}=h/(2e)$. In the last line of this equation, we have defined the vector $\vec{\phi }={\{{\phi }_{0},{\phi }_{1},\ldots ,{\phi }_{N},{\phi }_{{qb}}\}}^{T}$ of length $N+2$, and the capacitance and inductance matrices ${\boldsymbol{C}}$ and ${\boldsymbol{L}}$ such that equation (4) reproduces equation (3).

Ignoring the nonlinear term proportional to E_J for the moment, this Lagrangian leads to the Euler–Lagrange equation

$$\begin{eqnarray}&&\ddot{\vec{\phi }}=-{{\boldsymbol{C}}}^{-1}{\boldsymbol{L}}\vec{\phi }\equiv {{\boldsymbol{\Omega }}}^{2}\vec{\phi },\end{eqnarray} \tag{ 5 }$$

which has the qubit mode ${\vec{\phi }}_{{qb}}={\{0,\ldots ,0,{\phi }_{{qb}}\}}^{T}$ as one of the eigenvectors. Since the qubit's Josephson energy is not included in the matrix ${{\boldsymbol{\Omega }}}^{2}=-{{\boldsymbol{C}}}^{-1}{\boldsymbol{L}}$ of size $N+2\times N+2$, this mode has zero eigenvalue and can easily be identified. A convenient basis to treat the array and the qubit separately is obtained by finding the eigenvectors of the $N+1\times N+1$ block matrix of ${{\boldsymbol{\Omega }}}^{2}$ that does not relate to ${\vec{\phi }}_{{qb}}$. We refer to these eigenvectors as ${\vec{v}}_{k}$, such that the flux across the array is given as $\vec{\phi }(t)={\sum }_{k}{\phi }_{k}(t){\vec{v}}_{k}$ with the time-dependence written explicitly. With this approach, the array modes are already renormalized by the qubit capacitances, C_q and C_s. In contrast to the modes of a transmission line resonator, notice that the mode vectors, ${\overrightarrow{v}}_{k}$, are not orthogonal, i.e., ${\overrightarrow{v}}_{k}^{T}{\overrightarrow{v}}_{l}\ne 0$. These modes are, however, linearly independent such that each term of $\vec{\phi }(t)$ decouple from each other in the time evolution, equation (5), and thus constitute the normal modes of the array. In the basis of these eigenmodes, the Lagrangian in equation (4) takes the simple form

$$\begin{eqnarray}{ \mathcal L } & = & \displaystyle \sum _{k}\left[\displaystyle \frac{{C}_{k}}{2}{\dot{\phi }}_{k}^{2}-\displaystyle \frac{1}{2{L}_{k}}{\phi }_{k}^{2}-{C}_{q}{v}_{k}(0)\,{\dot{\phi }}_{k}{\dot{\phi }}_{{qb}}\right]\\ & & +\displaystyle \frac{{C}_{q}+{C}_{s}}{2}{\dot{\phi }}_{{qb}}^{2}+{E}_{J}\cos ({\phi }_{{qb}}\,/{\varphi }_{0}),\end{eqnarray} \tag{ 6 }$$

where ${v}_{k}(n)$ denotes the nth entry of the eigenvector ${\vec{v}}_{k}$. In the above expression, we have defined the mode capacitance C_k and mode inductance L_k as

$$\begin{eqnarray}&&{C}_{k}={\vec{v}}_{k}^{\,T}{\boldsymbol{C}}{\vec{v}}_{k},\qquad {L}_{k}^{-1}={\vec{v}}_{k}^{\,T}{\boldsymbol{L}}{\vec{v}}_{k}.\end{eqnarray} \tag{ 7 }$$

With these definitions, the eigenmode frequencies take the usual form ${\omega }_{k}=1/\sqrt{{L}_{k}{C}_{k}}$.

To obtain the associated Hamiltonian, we first identify the conjugate variables

$$\begin{eqnarray}&&{q}_{k}=\displaystyle \frac{\partial { \mathcal L }}{\partial {\dot{\phi }}_{k}}={C}_{k}{\dot{\phi }}_{k}-{C}_{q}{v}_{k}(0){\dot{\phi }}_{{qb}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{q}_{{qb}}=({C}_{q}+{C}_{s}){\dot{\phi }}_{{qb}}-\displaystyle \sum _{k}{C}_{q}{v}_{k}(0){\dot{\phi }}_{k}.\end{eqnarray} \tag{ 9 }$$

Introducing $\tilde{q}$ and $\tilde{\phi }$ as the row vectors of entries ${q}_{k({qb})}$ and ${\phi }_{k({qb})}$, the above expressions can be written in compact vector form as

$$\begin{eqnarray}&&\tilde{q}=\tilde{{\boldsymbol{C}}}\dot{\tilde{\phi }}.\end{eqnarray} \tag{ 10 }$$

We also define $\tilde{{\boldsymbol{L}}}$, the diagonal matrix of matrix elements $1/{L}_{k}$. Using this notation a Legendre transformation is performed and the Hamiltonian reads

$$\begin{eqnarray}&&H={\tilde{q}}^{T}\displaystyle \frac{{\tilde{{\boldsymbol{C}}}}^{-1}}{2}\tilde{q}+{\tilde{\phi }}^{T}\displaystyle \frac{\tilde{{\boldsymbol{L}}}}{2}\tilde{\phi }-{E}_{J}\cos {\phi }_{{qb}},\end{eqnarray} \tag{ 11 }$$

with the capacitances and inductances for mode k given by the diagonal entries of the matrices,

$$\begin{eqnarray}&&{\tilde{C}}_{k}^{-1}={\tilde{{\boldsymbol{C}}}}_{[k,k]}^{-1},\qquad {\tilde{L}}_{k}^{-1}={\tilde{{\boldsymbol{L}}}}_{[k,k]}.\end{eqnarray} \tag{ 12 }$$

The Hamiltonian of equation (11) can be expressed as the sum of a qubit Hamiltonian, H_qb, an array Hamiltonian, H_array, and their coupling, H_c. The qubit Hamiltonian takes the standard form

$$\begin{eqnarray}&&{H}_{{qb}}=4{E}_{C}{\hat{n}}^{2}-{E}_{J}\cos ({\phi }_{{qb}}/{\varphi }_{0}),\end{eqnarray} \tag{ 13 }$$

where $-2e\hat{n}={\hat{q}}_{{qb}}$ and ${E}_{C}={e}^{2}/(2{\tilde{C}}_{{qb}})$. In the transmon regime, ${E}_{J}/{E}_{c}\gg 1$, this Hamiltonian can be approximated as [23]

$$\begin{eqnarray}&&{H}_{{qb}}\approx {\omega }_{a}{b}^{\dagger }b-\displaystyle \frac{{E}_{C}}{2}{b}^{\dagger }{b}^{\dagger }{bb},\end{eqnarray} \tag{ 14 }$$

with the transmon frequency ${\omega }_{a}=\sqrt{8{E}_{J}{E}_{C}}-{E}_{C}$ and where we have introduced the operators

$$\begin{eqnarray}&&{\phi }_{{qb}}={\varphi }_{0}{\left(\displaystyle \frac{2{E}_{C}}{{E}_{J}}\right)}^{1/4}({b}^{\dagger }+b),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{n}_{{qb}}={\rm{i}}{\left(\displaystyle \frac{{E}_{J}}{32{E}_{C}}\right)}^{1/4}({b}^{\dagger }-b).\end{eqnarray} \tag{ 16 }$$

While the above form is useful in simplifying analytical expressions, all numerical calculations in this paper are based on the exact diagonalization of equation (13).

Expressing the mode operators q_k and ${\phi }_{k}$ in terms of the creation (annihilation) operators ${a}_{k}^{\dagger }$ (a_k) for mode k of the array

$$\begin{eqnarray}&&{q}_{k}={\rm{i}}\sqrt{\displaystyle \frac{{\tilde{\omega }}_{k}{\tilde{C}}_{k}}{2}}({a}_{k}^{\dagger }-{a}_{k}),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\phi }_{k}=\sqrt{\displaystyle \frac{{\tilde{\omega }}_{k}{\tilde{L}}_{k}}{2}}({a}_{k}^{\dagger }+{a}_{k}),\end{eqnarray} \tag{ 18 }$$

the array Hamiltonian takes the standard from ${H}_{{\rm{array}}}={\sum }_{k}{\tilde{\omega }}_{k}{a}_{k}^{\dagger }{a}_{k}$. In this expression, the mode frequencies are ${\tilde{\omega }}_{k}=1/({Z}_{k}{\tilde{C}}_{k})$ where ${Z}_{k}=\sqrt{{\tilde{L}}_{k}/{\tilde{C}}_{k}}$ is the characteristic impedance of mode k [27, 28]. The frequencies ${\tilde{\omega }}_{k}$ differ slightly from ${\omega }_{k}$ due to the off-diagonal elements of $\tilde{{\boldsymbol{C}}}$. As can be seen from the first term of equation (11), these terms also causes a small coupling between the array modes. This mode-mode coupling is due to the qubit-array interaction and is analogous to a multi-mode A²-term [31, 32]. Omitting array junctions nonlinearities, the Hamiltonian then takes the form

$$\begin{eqnarray}H & = & {H}_{{qb}}+{H}_{{\rm{array}}}+\displaystyle \sum _{k}{g}_{k}({b}^{\dagger }-b)({a}_{k}^{\dagger }-{a}_{k})\\ & & +\displaystyle \sum _{k\ne l}{G}_{{kl}}({a}_{k}^{\dagger }-{a}_{k})({a}_{l}^{\dagger }-{a}_{l}).\end{eqnarray} \tag{ 19 }$$

In the transmon regime ${E}_{J}/{E}_{C}\gg 1$, the qubit-array coupling strength takes the form

$$\begin{eqnarray}&&{g}_{k}={\left(\displaystyle \frac{8{E}_{J}}{{E}_{C}}\right)}^{\tfrac{1}{4}}\sqrt{\displaystyle \frac{1}{2{Z}_{k}}}e\ {\tilde{{\boldsymbol{C}}}}_{[k,{qb}]}^{-1}.\end{eqnarray} \tag{ 20 }$$

As expected, we find that ${g}_{k}/{\tilde{\omega }}_{k}\propto \sqrt{{Z}_{k}}$ [18, 27]. As already mentioned, in addition to a qubit-array coupling, equation (19) also contains a mode-mode interaction given by

$$\begin{eqnarray}&&{G}_{{kl}}=\sqrt{\displaystyle \frac{1}{2{Z}_{k}}}\sqrt{\displaystyle \frac{1}{2{Z}_{l}}}\ {\tilde{{\boldsymbol{C}}}}_{[k,l]}^{-1}/2.\end{eqnarray} \tag{ 21 }$$

In practice the frequency difference between modes is such that ${\tilde{\omega }}_{l}-{\tilde{\omega }}_{k}\gtrsim 100{G}_{{kl}}$, evaluated using the parameters used in section 3. Due to the small magnitude of these G_kl, we can neglect their renormalization of the mode-frequencies.

To finalize the derivation of the system Hamiltonian, we now include the array junction nonlinearities following the approach of [28, 29]. Taking advantage of the weak nonlinearity of these junctions, we consider only the fourth order expansion of the cosine potential of each array junction leading to the nonlinear potential

$$\begin{eqnarray}&&{U}_{{nl}}=-\displaystyle \frac{1}{24{L}_{J}{\phi }_{0}}\displaystyle \sum _{n=0}^{N-1}{({\phi }_{n}-{\phi }_{n+1})}^{4}.\end{eqnarray} \tag{ 22 }$$

Expressing this in the eigenmode basis, using the mode creation and annihilation operators, and dropping all rotating terms, this leads to the additional term in the Hamiltonian of equation (11) [28]

$$\begin{eqnarray}&&{H}_{K}=\displaystyle \sum _{{kl}}{K}_{{kl}}{a}_{k}^{\dagger }{a}_{k}{a}_{l}^{\dagger }{a}_{l},\end{eqnarray} \tag{ 23 }$$

where the self- (K_kk) and cross-Kerr (K_kl) coefficients can be expressed as

$$\begin{eqnarray}&&{K}_{{kl}}=-\displaystyle \frac{2-{\delta }_{{kl}}}{4{L}_{J}{\varphi }_{0}^{2}}\displaystyle \frac{{\tilde{\omega }}_{k}{\tilde{L}}_{k}}{2}\displaystyle \frac{{\tilde{\omega }}_{l}{\tilde{L}}_{l}}{2}\displaystyle \sum _{n=0}^{N-1}{\rm{\Delta }}{\phi }_{k}{(n)}^{2}{\rm{\Delta }}{\phi }_{l}{(n)}^{2},\end{eqnarray} \tag{ 24 }$$

with ${\rm{\Delta }}{\phi }_{k}(n)\equiv {v}_{k}(n)-{v}_{k}(n+1)$. It follows immediately that ${K}_{{kl}}=-2\sqrt{{K}_{{kk}}{K}_{{ll}}}$ [29].

We note that the Hamiltonian of equation (19) was obtained by first finding the qubit-renormalized array modes which were used as a convenient basis. With this approach, the modes already take into account the capacitances C_s and C_q which may be much larger than the array capacitances and hence significantly change the mode structure. With this choice, the mode-mode-couplings , G_kl, are then very small and can be ignored. Another approach to find the system Hamiltonian would be to first diagonalize the Lagrangian without coupling to ${\dot{\phi }}_{{qb}}$ and then reintroduce this coupling. Such an approach would lead to much larger mode-mode-coupling which would then have to be taken into account by exact diagonalization. Both approaches, in the end, lead to equivalent coupling strengths between the array modes and the qubit, g_k.

Before concluding this section, we note that the distinction between the resonator and the transmon in the system Hamiltonian may seem artificial. After all, the split into a transmon degree of freedom and array degrees of freedom is unnecessary to calculate the eigenfrequencies of the combined system. This distinction is, however, useful since one of the modes of the total system, the qubit mode, inherits the most from the transmon's large nonlinearity. This can be made more apparent by replacing the qubit junction by a SQUID. For symmetric junctions, this leads to the replacement ${E}_{J}\to {E}_{J}\cos ({{\rm{\Phi }}}_{x}/2{\varphi }_{0})$, with ${{\rm{\Phi }}}_{x}$ the external flux, in the qubit Hamiltonian H_qb. In this situation, the qubit mode is widely flux tunable while the array modes have, following our treatment, no explicit dependence on flux. This also affects the qubit-array coupling which, in the transmon regime, now takes the form

$$\begin{eqnarray}&&{g}_{k}({{\rm{\Phi }}}_{x})\approx \cos {({{\rm{\Phi }}}_{x}/2{\varphi }_{0})}^{1/4}\ {g}_{k}({{\rm{\Phi }}}_{x}=0),\end{eqnarray} \tag{ 25 }$$

with ${g}_{k}({{\rm{\Phi }}}_{x}=0)$ given by equation (20). As will be explored below, replacing the qubit junction by a SQUID also provides a tool for initiating dynamics in the system.

To investigate how strongly the transmon can be coupled to the array, we now focus on the lowest mode of the resonator, k = 0. Indeed, this mode is expected to have the largest zero-point fluctuations as characterized by ${g}_{0}/{\omega }_{0}\propto \sqrt{{Z}_{0}}$ (from now on we write ${\omega }_{k}$, however, still referring to ${\tilde{\omega }}_{k}$ calculated in section 2). To reach a large value of Z₀, the array junctions must be of large Josephson inductance L_J and of small capacitance to ground C₀. Moreover, since g₀ naturally depends on the coupling capacitor, C_q, it is also useful to make this capacitance large. However, a change in C₀ and C_q does not only change g₀, but it also influences other system parameters such as the transmon anharmonicity, E_C, the transmon frequency, ${\omega }_{a}$, the mode frequencies, ${\omega }_{k}$, and the other mode couplings, g_k. Our approach to maximize the coupling strength is thus to fix the qubit anharmonicity E_C and the mode frequency ${\omega }_{0}$ for fixed values of the capacitance C₀. The coupling g₀ is then optimized numerically by varying the rest of the system parameters. As will be clear below, C_q is not part of this optimization but we will varied to find an explicit dependence of g₀ on this capacitance. This approach does not guarantee the globally maximal coupling strength, but it is sufficient to identify parameters that yield a transmon in the ultrastrong coupling regime.

For the numerical examples presented below, we fix the first mode frequency to ${\omega }_{0}=2\pi \times 2\,\mathrm{GHz}$, while the transmon anharmonicity is fixed to ${E}_{C}=2\pi \times 300\,\mathrm{MHz}$ with Josephson energy ${E}_{J}/{E}_{C}=50$. For a resonator mode at this frequency cooled to 15 mK [24], the thermal population is as small as ∼ 0.0016 photons. On the other hand, larger effective temperatures have been reported in some experiments [33]. We note that, even at 50 mK, the thermal population is, however, only still of 0.14 photons. With a small value of this mode frequency, it is possible to reach a large ${g}_{0}/{\omega }_{0}$ ratio even with a moderate value of ${g}_{0}/2\pi \sim 1\,\mathrm{GHz}$. In turn, this means that non-adiabatic changes of parameters are possible with realistic flux modulations, allowing for the observation of ultrastrong dynamics. With these choices, figure 2 shows the results of a numerical optimization of the coupling strength as a function of C_q and for different values of C₀. As expected, increasing C_q leads to an increase of g₀. The observed oscillations in the coupling strength are due to local maxima in the numerical optimization. While the value of the coupling capacitor used here is larger than used in a previous realization of the ultrastrong coupling regime [34], it is comparable to capacitances used in transmon qubits [33]. More importantly, despite the large coupling capacitance, the parameters that we propose using here leads to the standard ${E}_{J}/{E}_{C}$ ratio of transmons. The results of figure 2 highlight that it is possible to reach the ultrastrong coupling regime with a large range of parameters. Finally, we note that the mode frequency ${\omega }_{0}$ was chosen to be small, but still large enough to avoid important thermal photon population.

Zoom In Zoom Out Reset image size

While generating non-trivial dynamics is the objective here, it also important to have a readout mechanism to probe this dynamic. Because of the photon-number dependent frequency shift resulting from cross-Kerr coupling, it is possible to use a second mode to probe the photon population of the fundamental mode [24]. Therefore, another design objective is to have a large cross-Kerr coupling between modes.

To reach these objectives, we take as parameters:

$$\begin{eqnarray*}N & = & 145,\qquad \ \ {L}_{J}=1.5\ \mathrm{nH},\\ {C}_{0} & = & 0.1\ \mathrm{fF},\qquad {C}_{q}=85\ \mathrm{fF},\\ {C}_{i} & = & 26\ \mathrm{fF},\qquad {C}_{J}=30\ \mathrm{fF},\\ {C}_{e} & = & 72\ \mathrm{fF},\qquad {C}_{s}=63\ \mathrm{fF},\end{eqnarray*}$$

which lead to

$$\begin{eqnarray*}{\omega }_{0}/2\pi & = & 2\,\ \mathrm{GHz},\quad \qquad \quad {g}_{0}/{\omega }_{0}=0.61,\\ {\omega }_{1}/2\pi & = & 8.8\,\ \mathrm{GHz},\,\,\,\qquad \,\ \ {g}_{1}/{\omega }_{1}=0.11,\\ {\omega }_{2}/2\pi & = & 14.45\ \ \mathrm{GHz},\,\quad \quad \,{g}_{2}/{\omega }_{2}=0.04,\\ {K}_{00}/2\pi & = & -0.03\,\ \mathrm{MHz},\qquad {K}_{22}/2\pi =-2.46\,\ \mathrm{MHz},\\ {K}_{02}/2\pi & = & -0.54\,\ \mathrm{MHz},\qquad {\omega }_{a}/2\pi =5.7\,\ \mathrm{GHz},\end{eqnarray*}$$

with ${\omega }_{a}$ the transmon qubit frequency at ${{\rm{\Phi }}}_{x}=0$. Furthermore we take the resonator decay rate $\kappa =2\pi \times 50\,\mathrm{kHz}$, corresponding to the losses observed in [24], which also includes losses induced by phase slips in the junction array. The qubit decay rate and the pure dephasing rate are taken as $\gamma ={\gamma }_{\phi }=2\pi \times 50\,\mathrm{kHz}$, values that are routinely observed for flux tunable transmons [35]. With these choices, the 0th mode is well within the ultrastrong coupling regime while the 1st mode is on the edge of that regime. Moreover, the 2nd mode is both outside the ultrastrong coupling regime and is far-detuned from the qubit. As a result, the dispersive coupling of that mode and any higher mode to the qubit is vanishingly small. On the other hand, as desired the 2nd mode has a significant cross-Kerr coupling to the 0th mode allowing for photon population readout.

Before moving to the dynamics of the system, it is instructive to also consider the spectral properties of the coupled transmon-array system by numerical diagonalizing of Hamiltonian equation (19). In figure 3(a), the eigenenergies for the transmon and the fundamental mode, k = 0, is shown as a function of the flux through the transmon SQUID, ${{\rm{\Phi }}}_{x}$, with the parameters listed above. The eigenmodes are highly hybridized due to the ultrastrong coupling and, as such, each eigenstate consists of both a substantial transmon and array component. To gain a more qualitative understanding of these results, we present in figures (b) and (c) a similar spectrum, but with the transmon replaced by a two-level system and a harmonic oscillator, respectively. The two-level system model of panel (b) is constructed by truncating the transmon Hilbert space to only two levels, while the harmonic oscillator is chosen such that its transition frequency is that of the transmon's 0-1 separation. We notice that, while the transmon hybridizes the energy levels in a similarly way as to the harmonic oscillator, the transmon spectrum includes a qubit-like avoided crossing indicated with the blue circles in figure 3. Moreover, the crossings highlighted by the red ellipse in harmonic oscillator spectrum are not present for the transmon. Therefore, while the transmon is not a qubit (i.e. a two-level system) due to its small anharmonicity, it pertains matter-like features similar to a two-level system.

Zoom In Zoom Out Reset image size

In this section, we present numerical results of the dynamics for the system with the above parameters and in the presence of damping. Because of the breakdown of the rotating-wave approximate in the ultrastrong coupling regime, it is not possible to use the standard quantum optics master equation [36]. We instead use a master equation derived in the instantaneous eigenbasis $\{| j(t)\rangle \}$ of the full system Hamiltonian including Kerr nonlinearity. Following [36], this master equation reads

$$\begin{eqnarray}\dot{\rho } & = & \displaystyle \frac{-{\rm{i}}}{{\hslash }}[H,\rho ]+\displaystyle \sum _{j,k\ne j}{{\rm{\Gamma }}}_{\phi }^{{jk}}{ \mathcal D }[| j\rangle \langle k| ]\rho \\ & & +\displaystyle \sum _{j,k\gt j}({{\rm{\Gamma }}}_{\kappa }^{{jk}}+{{\rm{\Gamma }}}_{\gamma }^{{jk}}){ \mathcal D }[| j\rangle \langle k| ]\rho \\ & & +{ \mathcal D }\left[\displaystyle \sum _{j}{{\rm{\Phi }}}^{j}| j\rangle \langle j| \right]\rho ,\end{eqnarray} \tag{ 26 }$$

with ${ \mathcal D }[O]\rho =O\rho {O}^{\dagger }-\tfrac{1}{2}({O}^{\dagger }O\rho +\rho {O}^{\dagger }O)$. This equation describes incoherent transitions and dephasing of the system eigenstates with the relaxation rates

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\gamma }^{{jk}}=\gamma \,| \langle j| {b}^{\dagger }+b| k\rangle {| }^{2},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\kappa }^{{jk}}=\kappa \,| \langle j| {a}_{0}^{\dagger }+{a}_{0}| k\rangle {| }^{2},\end{eqnarray} \tag{ 28 }$$

the dephasing rates

$$\begin{eqnarray}&&{{\rm{\Phi }}}^{j}=\sqrt{\displaystyle \frac{{\gamma }_{\phi }}{2}}| \langle j| {b}^{\dagger }b| j\rangle {| }^{2},\end{eqnarray} \tag{ 29 }$$

and the dephasing-induced relaxation rates

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\phi }^{{jk}}=\displaystyle \frac{{\gamma }_{\phi }}{2}| \langle j| {b}^{\dagger }b| k\rangle {| }^{2}.\end{eqnarray} \tag{ 30 }$$

In the above expressions, a₀ (${a}_{0}^{\dagger }$) refers to the fundamental mode annihilation (creation) operator and b (${b}^{\dagger }$) to the qubit lowering (raising) operator. With these forms for the rates, the equilibrium state of equation (26) is the ground state of the coupled system [36]. In contrast, the quantum optics master equation would bring the system to the ground state of the uncoupled system, a state which is far from the true ground state in the ultrastrong coupling regime.

4.1. Non-adiabatic generation of photons

As already mentioned, an important feature of the Jaynes–Cummings Hamiltonian is that its ground state is that of the uncoupled system. As a result, the nature of this ground state does not change with system parameters. In other words, if prepared in its ground state, a system described by the Jaynes–Cummings Hamiltonian will remain in the vacuum state under parametric modulations.

In contrast, the ground state, $| j=0\rangle$, of the Rabi Hamiltonian can be approximated as [36]

$$\begin{eqnarray}&&| j=0\rangle \approx \left(1-\displaystyle \frac{{{\rm{\Lambda }}}^{2}}{2}\right)| 00\rangle -{\rm{\Lambda }}| 11\rangle +\xi \sqrt{2}| 02\rangle \end{eqnarray} \tag{ 31 }$$

to second order in ${\rm{\Lambda }}=g/({\omega }_{a}+{\omega }_{r})$ and with $\xi =g{\rm{\Lambda }}/2{\omega }_{r}$. On the right-hand-side of this expression, the first index in the states refers to the qubit and the second to the photon number. Equation (31) makes it clear that the ground state of H_Rabi depends on the system parameters and, moreover, has a finite average photon number. Since the master equation equation (26) relaxes the system back to $| j=0\rangle$, these photons do not decay out of the cavity and are consequently difficult to observe.

Here we propose to take advantage of the dependence of $| j=0\rangle$ on the system parameters to observe a signature of these photons. Indeed, a non-adiabatic change of the system parameters should lead to a change of the average photon population under H_Rabi while it should have no effect under H_JC. As alluded to earlier, this photon population can then be probed by taking advantage of the cross-Kerr coupling between the array modes. For the photon population to change under parametric modulations, this modulation must, however, be non-adiabatic. This is possible in this system and with the parameters of section 3 because of the small mode frequency ${\omega }_{0}$ and therefore the reasonably small g₀ required to reach ultrastrong coupling.

To realize this, we modulate the flux through the transmon's SQUID loop as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{x}(t)={{\rm{\Phi }}}_{x}^{a}\cos ({\omega }_{d}t)\end{eqnarray} \tag{ 32 }$$

to induce non-adiabatic dynamics [37]. To reach measurable photon populations, large flux modulations $\gtrsim 0.1{{\rm{\Phi }}}_{0}$ are required. While this modulation amplitude is larger than what is typically used in flux-pumped Josephson parametric amplifiers, similar amplitudes have already been demonstrated experimentally [38].

Because of the change in system parameters under this flux modulation, the overlap between the instantaneous ground state at a given time, $| j=0(t)\rangle$, and the $j^{\prime}$th excited state at a later time $t^{\prime}$, $| j^{\prime} (t^{\prime} )\rangle$ will in general be non-zero,

$$\begin{eqnarray}&&\langle j=0(t)| j^{\prime} (t^{\prime} )\rangle \ne 0,\end{eqnarray} \tag{ 33 }$$

a result that holds only when the RWA is not valid. This implies that flux modulations can excite the system away from the ground state. An example of this non-adiabatic dynamics is presented in figure 4(a) which shows the photon population as a function of time as obtained by numerical integration of equation (26) with the parameters of section 3 and a modulation frequency of ${\omega }_{d}=2\pi \times 1.5\,\mathrm{GHz}$. In these simulations, the system was first initialized in the ground state $| j=0\rangle$. Importantly, the drive frequency does not correspond to a resonance frequency of the coupled system and is therefore not expected to directly drive specific system transitions. Despite this, a consequential photon population is observed for ${g}_{0}/{\omega }_{0}=0.61$ (dark-blue line). On the other hand, a weaker coupling of ${g}_{0}/{\omega }_{0}=0.1$ (light blue line) for which the Jaynes–Cummings Hamiltonian is expected to be a good approximation shows a much smaller average photon population.

Zoom In Zoom Out Reset image size

To further illustrate this point, figure 5 shows the time-averaged photon number as a function of the modulation frequency ${\omega }_{d}$ and for different drive strengths. Again, we observe that, for the small coupling strengths, very few photons are generated and this occurs only at well-defined resonances. In the ultrastrong coupling regime photons are, however, generated for a large range of frequencies. For the strongest coupling of ${g}_{0}/{\omega }_{0}=0.61$ (dark blue line), photons are observed for all drive frequencies. Moreover, for some specific frequencies, we observe an enhancement of the photon generation. Due to the complexity of the transmon-array spectrum, these features can be difficult to identify as specific transitions in the spectrum. Additionally, we notice that the generated photon number is significantly larger than the expected thermal population of the array. The fact that any modulation frequency leads to photons in the array confirms that photons are not generated by directly exciting a transition of the static system, but are due to the non-adiabatic change of the ground state. An analogy can be drawn to multi-passage Landau-Zener transitions [39]. These transitions appear when the parameters of a two-level system is changed in a non-adiabatic fashion through an avoided crossing. A similar effect is observed here with a non-adiabatic change in the ground state. The transmon-array system has, however, a complex level structure where the Landau-Zener results cannot be explicitly applied. However, as illustrated in figure 3, avoided crossings appear both in the spectrum obtained for a transmon and a two-level system. In both cases this helps towards the generation of photons in the array. In figure 4(b), we see that in the two-level system plus array, the photon number admits a much richer dynamics due to the stronger nonlinearity and broader avoided crossing. In figure 4(b), we also compare the situation to a purely harmonic system. As illustrated in figure 3(c), this system does not shown avoided crossings and the modulation can here be considered as a simple change of boundary condition of the array mode. This modulation is then comparable to the dynamical Casimir effect [38]. In practice, we see that the effect is much weaker than the non-adiabatic generation of photons that are observed for the ultrastrong light–matter interaction with a transmon and a two-level system. Finally, we emphasize the qubit-like behavior of the transmon by pointing out the populations of the transmon, shown in figure 4(c). We observe the small population beyond the first two level, with less than 0.1% above the first 4 levels. This emphasize the few-level matter-like behavior of the transmon.

Zoom In Zoom Out Reset image size

4.2. Photon population measurement

To measure the photon population in the ultrastrongly coupled mode k = 0, we take advantage of the cross-Kerr coupling between modes k = 0 and 2. This coupling was already used experimentally to characterize a junction array [24]. Ignoring the other array modes, this coupling takes the form

$$\begin{eqnarray}&&{H}_{K}={K}_{02}{a}_{0}^{\dagger }{a}_{0}\,{a}_{2}^{\dagger }{a}_{2},\end{eqnarray} \tag{ 34 }$$

with ${K}_{02}=-4\sqrt{{K}_{00}{K}_{22}}$. Photon population in mode 0 will shift the second mode frequency by ${K}_{02}{a}_{0}^{\dagger }{a}_{0}$, a shift that can be resolved by probing mode 2.

The general approach is now to apply a coherent drive, ${H}_{p}={\varepsilon }_{p}({a}_{2}+{a}_{2}^{\dagger })$, on resonance with the probe mode via the input port C_i (see figure 1). The signal reflected from this port is then continuously monitored. Similarly to dispersive qubit readout [40], the photon number, $\langle {a}_{0}^{\dagger }{a}_{0}\rangle$, can be determined by homodyne measurement of the field amplitude a₂. The integrated homodyne signal can be expressed as

$$\begin{eqnarray}&&{M}_{{\varepsilon }_{p}}=\sqrt{{\kappa }_{2}}{\int }_{\tau }^{\tau +{T}_{m}}[{a}_{\mathrm{out}}(t)+{a}_{\mathrm{out}}^{\dagger }(t)]\,{\rm{d}}{t}\end{eqnarray} \tag{ 35 }$$

with T_m the integration time and τ the initial time of the integration. In this expression, ${a}_{\mathrm{out}}(t)=\sqrt{{\kappa }_{2}}{a}_{2}(t)+{a}_{\mathrm{in}}(t)$ is the output field [41], with a_in the input noise of the vacuum respecting $[{a}_{\mathrm{in}}(t),{a}_{\mathrm{in}}^{\dagger }(t^{\prime} )]=\delta (t^{\prime} -t)$, and ${\kappa }_{2}$ the decay rate of mode k = 2. The ability for such a measurement to distinguish the state from the system with no flux modulation, ie. no ultrastrong dynamics, is captured by the signal-to-noise ratio (SNR). Following [42, 43], the SNR can be expressed as

$$\begin{eqnarray}&&{\mathrm{SNR}}_{{T}_{m}}=\displaystyle \frac{| \langle {M}_{{\varepsilon }_{p}}\rangle -\langle {M}_{0}\rangle | }{\sqrt{\langle {\tilde{M}}_{{\varepsilon }_{p}}^{2}\rangle +\langle {\tilde{M}}_{0}^{2}\rangle }},\end{eqnarray} \tag{ 36 }$$

with ${\tilde{M}}_{{\varepsilon }_{p}}={M}_{{\varepsilon }_{p}}-\langle {M}_{{\varepsilon }_{p}}\rangle$ and M₀ corresponds to the same measurement without flux-modulation, ${{\rm{\Phi }}}_{x}=0$.

As illustrated in figure 6, two approaches are considered. In the first approach, depicted in panel (a), the probe field is monitored while continuously modulating the qubit flux. For simplicity, the nonlinearity of the probe mode is ignored and the photon number ${a}_{0}^{\dagger }{a}_{0}$ is taken to be a classical number. Then, the equation of motion for a₂ reads

$$\begin{eqnarray}&&{\dot{a}}_{2}=-{{\rm{i}}{K}}_{02}\langle {a}_{0}^{\dagger }{a}_{0}\rangle {a}_{2}-\displaystyle \frac{{\kappa }_{2}}{2}{a}_{2}-{\rm{i}}{\varepsilon }_{p}+\sqrt{{\kappa }_{2}}{a}_{\mathrm{in}}\end{eqnarray} \tag{ 37 }$$

which as the steady state solution

$$\begin{eqnarray}&&{a}_{2}^{{\rm{s}}}=\displaystyle \frac{\sqrt{{\kappa }_{2}}{a}_{\mathrm{in}}-{\rm{i}}{\varepsilon }_{p}}{{{\rm{i}}{K}}_{02}\langle {a}_{0}^{\dagger }{a}_{0}\rangle +{\kappa }_{2}/2}.\end{eqnarray} \tag{ 38 }$$

To obtain a simple estimate for the SNR, we use the values of $\langle {a}_{0}^{\dagger }{a}_{0}\rangle$ oscillating between 0.2 and 0.8 shown in figure 4 and integrate the signal taking $\tau =100$ ns to go beyond the initial ring up dynamics. For the parameters presented in section 3, together with ${\varepsilon }_{p}=2\pi \times 2$ MHz and ${\kappa }_{2}=2\pi \times 0.35\,\mathrm{MHz}$, this yields a SNR larger than 1 for an integration time ${T}_{m}\approx {\kappa }_{2}^{-1}$. A larger SNR can be obtained by longer integration times, however, the ultrastrong dynamics will eventually dephase due the dephasing rates ${{\rm{\Phi }}}^{j}$. Using equation (38) the value of $\langle {a}_{0}^{\dagger }{a}_{0}\rangle$ is estimated and we recover, as desired, the numerical time-averaged results shown in figure 5. In this analysis we neglected the self-Kerr nonlinearity, ${K}_{22}=2\pi \times -2.4\,\mathrm{MHz}$, but in general similar results for the cross-Kerr probing can be obtained by including the nonlinearity in the analysis [44].

Zoom In Zoom Out Reset image size

An alternative method to map the dynamics shown in figure 4 is sketched in figure 6(b). In this approach, the qubit flux is modulated for a time τ around ${{\rm{\Phi }}}_{x}=0$ with an amplitude of ${{\rm{\Phi }}}_{x}=0.35{{\rm{\Phi }}}_{0}$. After this initial period of ultrastrong dynamics, the flux is rapidly increased to ${{\rm{\Phi }}}_{x}={{\rm{\Phi }}}_{0}/2$ in a time span of one full period of oscillation, $2\pi /{\omega }_{d}\approx 0.66$ ns. At that point, the qubit has a vanishingly small transition frequency and is uncoupled from the array, see equation (25). Now, with the coupling to the qubit absent, the population of mode a₀ simply decays to the vacuum state at a rate ${\kappa }_{0}$. Again, it is worth emphasizing that due to the choice of a small resonator frequency and, thus, low coupling, this change in flux is fast enough to maintain the photon number in the a₀ mode. Using the same parameters as in figure 4, we numerically integrate equation (36) and find, taking into account array damping, a maximal SNR${}_{{T}_{m}}\approx 0.5$ for a measurement time ${T}_{m}\approx 6$ μs. For larger measurement times, the signal will be dominated by noise because the photon population of the a₀ mode have decayed. This estimate is obtained from numerical integration including cross-Kerr coupling given by equation (34) and self-Kerr nonlinearities for both modes. As above, from the measured signal, the detuning of the probe mode a₂ from its bare frequency ${\omega }_{2}$ can now be inferred. Using the inferred probe detuning, the photon population in mode k = 0 can then be estimated. Therefore, by treating the average photon number as an unknown parameter, standard parameter estimation techniques [45, 46] can be used and the dynamics generated during the ultrastrong coupling is observed. We expect that a measurement of the dynamics can be obtained by only a few experimental runs for each value of τ.

In conclusion, we have shown that it is possible to reach the ultrastrong coupling regime of light–matter interaction by coupling a transmon qubit to a high impedance mode realized by an array of weakly nonlinear Josephson junctions. Using realistic system parameters, we find coupling strengths as large as 0.6 times the system frequency. By working with a low frequency mode of the array, this ultrastrong coupling is obtained for moderate values of the coupling. This is an important advantage of our proposal. Indeed, with this choice, we have shown that realistic modulations of the transmon parameters are sufficient to result in dynamics of the system that is distinctive of the ultrastrong coupling regime. Moreover, we have shown how this dynamic and the corresponding photon population can be probed by taking advantage of the multi-mode structure of the array. These results show the possibility to probe the complex dynamics of the ultrastrong coupling regime, opening a new window on this unconventional regime of quantum optics.

Although, these results were obtained using a transmon, the general idea applies to different types of artificial atoms. For example, it could be interesting to take advantage of the large anharmonicity of flux qubits [47] or of the rich level structure of the fluxonium [48] in probing the ultra-strong dynamics.

The authors are grateful to J Bourassa, N Didier, A Kamal and N Roch for discussions. Financial support from the Villum Foundation Center of Excellence, QUSCOPE, from the Danish Ministry of Higher Education and Science and by NSERC is acknowledged. This research was undertaken thanks in part to funding from the Canada First Research Excellence Fund.