Ultrastrong coupling dynamics with a transmon qubit

The interaction of light and matter is often described by the exchange of single excitations. When the coupling strength is a significant fraction of the system frequencies, the number of excitations are no longer preserved and that simple picture breaks down. This regime is known as the ultrastrong coupling regime and is characterized by non-trivial light-matter eigenstates and complex dynamics. In this work, we propose to use a an array Josephson junctions to increase the impedance of the light mode enabling ultrastrong coupling to a transmon qubit. We show that the resulting dynamics can be generated and probed by taking advantage of the multi-mode structure of the junction array. This proposal relies on the frequency tunability of the transmon and, crucially, on the use of a low frequency mode of the array, which allows for non-adiabatic changes of the ground state.


I. INTRODUCTION
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 highquality 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 ( = 1), where ω r is the cavity frequency, ω a the atomic frequency and g the electric-dipole coupling. In this expression, a (a † ) is the photon annihilation (creation) operator and σ 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 in 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 † + a), and the atomic dipole moment, d 0 σ x , where g = d 0 E 0 . The Jaynes-Cummings Hamiltonian is a good approximation to H Rabi when the coupling, g, is * Email: ctc@phys.au.dk smaller than the system frequencies, g ω a , ω r . In this situation, the fast rotating term a † σ + + aσ − appearing in the Rabi Hamiltonian can safely be dropped using the rotating wave approximation (RWA) and we recover Eq. (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][7][8][9][10][11][12][13]. Most significantly, while the ground state of H 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 nonzero average photon number. From a practical point of view, this regime could also be useful in the context of quantum information processing [14][15][16][17].
Superconducting quantum circuits form a promising platform to realize and study this novel light-matter coupling regime. In particular, 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][20][21][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][19][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 fo- Junction array Circuit representation of the system. An array of N Josphson junctions treated as a series of inductors with inductance LJ and capacitance CJ . These junctions also have a parasitic capacitance to ground, C0. The array is interaction via Cq with a transmon qubit characterized by the capacitance Cs and the Josephson energy EJ . The flux node at the transmon is denoted φ qb and the nodes of the array goes from φ0 to φN .
cus 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][25][26]. Together with its capacitance to ground, this superinductance plays the role of a multimode cavity. With g/ω r ∝ √ Z 0 , where Z 0 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. This dynamic 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 Sec. II the Hamiltonian of the system, taking into account the multi-mode structure of the array. In Sec. III, we identify parameters to reach the ultrastrong coupling regime. In Sec. IV the dynamics of the ultrastrong coupling regime is investigated. Finally, Sec. V concludes the paper.

II. TRANSMON COUPLED TO A JOSEPHSON JUNCTION ARRAY
We consider the circuit of Fig. 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 charac-terized by their Josephson inductance L J and junction capacitance C J . Following Refs. [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 0 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].
Following the standard approach [30], the circuit Lagrangian reads where φ n is the node flux of the nth island of the array, φ qb the node flux of the transmon's island and ϕ 0 = Φ 0 /2π with the magnetic flux quantum Φ 0 = h/(2e).
In the last line of this equation, we have defined the vector φ = {φ 0 , φ 1 , . . . , φ N , φ qb } T of length N + 2, and the capacitance and inductance matrices C and L such that Eq. (4) reproduces Eq. (3). Ignoring the non-linear term proportional to E J for the moment, this Lagrangian leads to the Euler-Lagrange equation¨ which has the qubit mode φ qb = {0, . . . , 0, φ qb } T as one of the eigenvectors. Since the qubit's Josephson energy is not included in the matrix Ω 2 = −C −1 L of size N +2×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×N +1 block matrix of Ω 2 that does not relate to φ qb . We refer to these eigenvectors as v k , such that the flux across the array is given as φ(t) = k φ k (t) 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 the basis of these eigenmodes, the Lagrangian Eq. (4) takes the simple form where v k (n) denotes the nth entry of the eigenvector v k . In the above expression, we have defined the mode capacitance C k and mode inductance L k as With these definitions, the eigenmode frequencies take the usual form ω k = 1/ √ L k C k . To obtain the associated Hamiltonian, we first identify the conjugate variables Introducingq andφ as the row vectors of entries q k (qb) and φ k (qb) , the above expressions can be written in compact vector form asq =Cφ.
We also defineL, the diagonal matrix of matrix elements 1/L k . Using this notation a Legendre transformation is performed and the Hamiltonian reads with the capacitances and inductances for mode k given by the diagonal entries of the matrices, The Hamiltonian of Eq. (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 where −2en =q qb and E C = e 2 /(2C qb ). In the transmon regime, E J /E c 1, this Hamiltonian can be approximated as [23] with the transmon frequency ω a = √ While the above form is useful in simplifying analytical expressions, all numerical calculations in this paper are based on the exact diagonalization of Eq. (13).
Expressing the mode operators q k and φ k in terms of the creation (annihilation) operators a † k (a k ) for mode k of the array the array Hamiltonian takes the standard from H array = kω k a † k a k . In this expression, the mode frequencies arẽ [27,28]. The frequenciesω k differ slightly from ω k due to the off-diagonal elements ofC. As can be seen from the first term of Eq. (11), these terms also causes a small coupling between the array modes. This mode-mode coupling is due to the qubitarray interaction and is analogous to a multi-mode A 2term [31,32]. Omitting array junctions nonlinearities, the Hamiltonian then takes the form In the transmon regime E J /E C 1, the qubit-array coupling strength takes the form As expected, we find that g k /ω k ∝ √ Z k [18,27]. As already mentioned, in addition to a qubit-array coupling, Eq. (19) also contains a mode-mode interaction given by In practice the frequency difference between modes is such thatω l −ω k 100G kl , evaluated using the parameters used in Sec. III. Due to the small magnitude of these G kl , we can neglect their renormalization of the modefrequencies.
To finalize the derivation of the system Hamiltonian, we now include the array junction nonlinearities following the approach of Refs. [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 non-linear potential 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 Eq. (11) [28] where the self-(K kk ) and cross-Kerr (K kl ) coefficients can be expressed as [29]. We note that the Hamiltonian of Eq. (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φ 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, leads 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 → E J cos(Φ x /2ϕ 0 ), with Φ 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 with g k (Φ x = 0) given by Eq. (20). As will be explored below, replacing the qubit junction by a SQUID also provides a tool for initiating dynamics in the system.

III. ULTRASTRONG COUPLING WITH A TRANSMON
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 /ω 0 ∝ √ Z 0 (from now on we write ω k , however, still referring toω k calculated in Sec. II). To reach a large value of Z 0 , the array junctions must be of large Josephson inductance L J and of small capacitance to ground C 0 . Moreover, since g 0 naturally depends on the coupling capacitor, C q , it is also useful to make this capacitance large. However, a change in C 0 and C q does not only change g 0 , but it also influences other system parameters such as the transmon anharmonicity, E C , the transmon frequency, ω a , the mode frequencies, ω 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 ω 0 for fixed values of the capacitance C 0 . The coupling g 0 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 0 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 ω 0 = 2π × 2 GHz, while the transmon anharmonicity is fixed to E C = 2π × 300 MHz with Josephson energy E J /E C = 50. With a small value of this mode frequency, it is possible to reach a large g 0 /ω 0 ratio even with a moderate value of g 0 /2π ∼ 1 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, Fig. 2 shows the results of a numerical optimization of the coupling strength as a function of C q and for different values of C 0 . As expected, increasing C q leads to an increase of g 0 . The observed oscillations in the coupling strength are due to local maxima in the numerical optimization. The results of Fig. 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 ω 0 was chosen to be small, but still large enough to avoid important thermal photon population.
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: N = 145, L J = 1.5 nH, C 0 = 0.1 fF, C q = 85 fF, C i = 26 fF, C J = 30 fF, C e = 72 fF, C s = 63 fF, which lead to with ω a the transmon qubit frequency at Φ x = 0. Furthermore we take the resonator decay rate κ = 2π×50 kHz, corresponding to the losses observed in Ref. [24]. The qubit decay rate and the pure dephasing rate are taken as γ = γ φ = 2π×50 kHz, values that are routinely observed for flux tunable transmons [33]. 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 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.

IV. DYNAMICS IN THE ULTRASTRONG COUPLING REGIME
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 [34]. We instead use a master equation derived in the instantaneous eigenbasis {|j(t) } of the full system Hamiltonian including Kerr nonlinearity. Following Ref. [34], this master equation readṡ with . This equation describes incoherent transitions and dephasing of the system eigenstates with the relaxation rates the dephasing rates and the dephasing-induced relaxation rates In the above expressions, a 0 (a † 0 ) refers to the fundamental mode annihilation (creation) operator and b (b † ) to the qubit lowering (raising) operator. With these forms for the rates, the equilibrium state of Eq. (26) is the ground state of the coupled system [34]. 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.

A. 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 , of the Rabi Hamiltonian can be approximated as [34] to second order in Λ = g/(ω a + ω r ) and with ξ = gΛ/2ω 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 Eq. (26) relaxes the system back to |j = 0 , 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 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 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 Sec. III because of the small mode frequency ω 0 and therefore the reasonably small g 0 required to reach ultrastrong coupling.
To realize this, we modulate the flux through the transmon's SQUID loop as to induce non-adiabatic dynamics [35]. To reach measurable photon populations, large flux modulations 0.1Φ 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 [36]. 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) , and the j th excited state at a later time t , |j (t ) will in general be non-zero, 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 dynamic is presented in Fig. 3 which shows the photon population as a function of time as obtained by numerical integration of Eq. (26) with the parameters of Sec. III and a modulation frequency of ω d = 2π × 1.5 GHz. In these simulations, the system was first initializing in the ground state |j = 0 . 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 /ω 0 = 0.61 (dark-blue line). On the other hand, a weaker coupling of g 0 /ω 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.
To further illustrate this point, Fig. 4 shows the timeaveraged photon number as a function of the modulation frequency ω 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 welldefined resonances. In the ultrastrong coupling regime photons are, however, generated for a large range of frequencies. For the strongest coupling of g 0 /ω 0 = 0.61 (dark blue line), photons are observed for all drive frequencies. This 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-  Fig. 3.
Cross-Kerr readout schemes to probe the photons generated by the ultrastrong coupling dynamics. In (a) we apply the modulation of the transmon continuously while probing a higher mode of the Kerr resonator. In (b) we only modulate for a time τ followed by a probing of a higher mode.
Zener transitions [37]. These transitions appear when the parameters of a two-level system is changed in a nonadiabatic 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.

B. 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 with K 02 = −4 √ K 00 K 22 . Photon population in mode 0 will shift the second mode frequency by K 02 a † 0 a 0 , a shift that can be resolved by probing mode 2.
The general approach is now to apply a coherent drive, H p = ε p (a 2 + a † 2 ), on resonance with the probe mode via the input port C i (see Fig. 1). The signal reflected from this port is then continuously monitored. Similarly to dispersive qubit readout [38], the photon number, a † 0 a 0 , can be determined by homodyne measurement of the field amplitude a 2 . The integrated homodyne signal can be expressed as with T m the integration time and τ the initial time of the integration. In this expression, a out (t) = √ κ 2 a 2 (t)+a in (t) is the output field [39], with a in the input noise of the vacuum respecting [a in (t), a † in (t )] = δ(t − t), and κ 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 Ref. [40,41], the SNR can be expressed as withM εp = M εp − M εp and M 0 corresponds to the same measurement without flux-modulation, Φ x = 0. As illustrated in Fig. 5, 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 a 0 is taken to be a classical number. Then, the equation of motion for a 2 readṡ which as the steady state solution To obtain a simple estimate for the SNR, we use the values of a † 0 a 0 oscillating between 0.2 and 0.8 shown in Fig. 3 and integrate the signal taking τ = 100 ns to go beyond the initial ring up dynamics. For the parameters presented in Sec. III, together with ε p = 2π × 2 MHz and κ 2 = 2π × 0.35 MHz, this yields a SNR larger than 1 for an integration time T m ≈ κ −1 2 . A larger SNR can be obtained by longer integration times, however, the ultrastrong dynamics will eventually dephase due the dephasing rates Φ j . Using Eq. (38) the value of a † 0 a 0 is estimated and we recover, as desired, the numerical time-averaged results shown in Fig. 4. In this analysis we neglected the self-Kerr nonlinearity, K 22 = 2π × −2.4 MHz, but in general similar results for the cross-Kerr probing can be obtained by including the nonlinearity in the analysis [42].
An alternative method to map the dynamics shown in Fig. 3 is sketched in Fig. 5(b). In this approach, the qubit flux is modulated for a time τ around Φ x = 0 with an amplitude of Φ x = 0.35Φ 0 . After this initial period of ultrastrong dynamics, the flux is rapidly increased to Φ x = Φ 0 /2 in a time span of one full period of oscillation, 2π/ω d ≈ 0.66 ns. At that point, the qubit has a vanishingly small transition frequency and is uncoupled from the array, see Eq. (25). Now, with the coupling to the qubit absent, the population of mode a 0 simply decays to the vacuum state at a rate κ 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 0 mode. Using the same parameters as in Fig. 3, we numerically integrate Eq. (36) and find, taking into account array damping, a maximal SNR Tm ≈ 0.5 for a measurement time T m ≈ 6 µs. For larger measurement times, the signal will be dominated by noise because the photon population of the a 0 mode have decayed. This estimate is obtained from numerical integration including cross-Kerr coupling given by Eq. (34) and self-Kerr nonlinearities for both modes. As above, from the measured signal, the detuning of the probe mode a 2 from its bare frequency ω 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 [43,44] 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 τ .

V. CONCLUSION
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 qubit parameters are sufficient to results 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.