---
abstract: 'We present Rabi oscillation measurements of a Nb/AlO$_\text{x}$/Nb dc superconducting quantum interference device (SQUID) phase qubit with a 100 ${\ensuremath{\mu {\ensuremath{\mathrm{m}}}}}^2$ area junction acquired over a range of microwave drive power and frequency detuning. Given the slightly anharmonic level structure of the device, several excited states play an important role in the qubit dynamics, particularly at high power. To investigate the effects of these levels, multiphoton Rabi oscillations were monitored by measuring the tunneling escape rate of the device to the voltage state, which is particularly sensitive to excited state population. We compare the observed oscillation frequencies with a simplified model constructed from the full phase qubit Hamiltonian and also compare time-dependent escape rate measurements with a more complete density-matrix simulation. Good quantitative agreement is found between the data and simulations, allowing us to identify a shift in resonance (analogous to the ac Stark effect), a suppression of the Rabi frequency, and leakage to the higher excited states.'
author:
- 'S. K. Dutta'
- 'Frederick W. Strauch'
- 'R. M. Lewis'
- Kaushik Mitra
- Hanhee Paik
- 'T. A. Palomaki'
- Eite Tiesinga
- 'J. R. Anderson'
- 'Alex J. Dragt'
- 'C. J. Lobb'
- 'F. C. Wellstood'
title: Multilevel effects in the Rabi oscillations of a Josephson phase qubit
---

Introduction {#SIntro}
============

Over the past decade, a variety of qubits based on superconducting Josephson junctions has been proposed and experimentally realized.[@Makhlin01a; @Devoret04a; @You05a; @Wendin07a; @Zagoskin07a] Among them is the phase qubit, which can take the form of a single current-biased junction,[@Ramos01a; @Yu02a] a flux-biased rf superconducting quantum interference device (SQUID),[@Simmonds04a] a low inductance dc SQUID,[@Claudon04a] a large inductance dc SQUID where one of the junctions is thought of as the qubit,[@Martinis02a] or a vortex in an annular junction.[@Wallraff03c] The behaviors of the devices are similar, as they all employ a large area junction, where the dynamics of the quantum mechanical phase difference across the junction (or the orientation of the vortex in the last case) are determined by a tilted washboard-like potential. The two lowest states in a well of this potential serve as the qubit basis. Among the possible benefits of the phase qubit are its relative insensitivity to charge and flux noise and its ability to operate over a wide range of parameters. Recent demonstrations include fast readout,[@Claudon04a] Rabi oscillations,[@Yu02a; @Martinis02a; @Claudon04a; @Strauch06a; @Lisenfeld07a] and simple capacitive coupling of two qubits as measured through spectroscopy,[@Xu05a] simultaneous state measurement,[@McDermott05a] and state tomography.[@Steffen06b]

In the phase qubit, higher excited states can impact the dynamics strongly, due to the nearly harmonic level structure within the potential well. The presence of these quantized energy levels has been clearly detected by monitoring the decay of a thermal population[@Silvestrini97a] and by microwave spectroscopy with single[@Martinis85a] and multiphoton[@Wallraff03a; @Strauch06a] transitions. Multilevel systems can, for example, be exploited for state initialization[@Valenzuela06a] and readout,[@Martinis02a] quantum logic gates[@Amin03a; @Strauch03a; @Wei08a] and algorithms,[@Ahn00a] and cryptography.[@Groblacher06a] However, when controlling the state of the phase qubit with a microwave current, off-resonant excitation of higher levels leads to leakage out of the desired qubit space and therefore a loss of quantum information.

Rabi oscillations serve as both a standard demonstration of quantum state manipulation and a diagnostic for decoherence and control fidelity. While these oscillations take on a simple form in a two-level system, several effects can occur in a multilevel system: the oscillations can distort, the higher states can become populated, and the basic resonance properties of the oscillations can undergo subtle shifts. Ultimately, to achieve fast and accurate control of the qubit state, all of these effects need to be carefully characterized. Progress in this area has been predominantly theoretical in nature, pointing out the conditions under which errors are introduced and methods to minimize their effects.[@Goorden03a; @Steffen03a; @Meier05a; @Amin06a; @Strauch06a; @Shevchenko07a] There has, however, been less direct experimental evidence confirming the validity of models upon which these predictions are based,[@Claudon04a; @Claudon07a; @Strauch06a; @Lucero08a] particularly with regard to leakage.

In this paper, we investigate the influence of the higher states of a phase qubit by examining the behavior of Rabi oscillations at high power, where their impact is greatest. In Sec. \[SPhase\], we describe the design of our qubit and a Hamiltonian that approximates its dynamics. Also discussed are various experimental details and the readout scheme that allows measurement of very small excited state populations. Section \[SStrong\] describes Rabi oscillations taken at a range of power and detuning from resonance, along with comparisons to a simple theory, while in Sec. \[SDecoherence\] we develop a more complete density-matrix model including the effects of decoherence and noise. Finally, Sec. \[SSummary\] contains a summary of the key points of our findings.

 SQUID Phase Qubit {#SPhase}
==================

![\[FDevice\]The dc SQUID phase qubit. (a) The qubit junction $J1$ (with critical current $I_{01}$ and capacitance $C_1$) is isolated from the current bias leads by an auxiliary junction $J2$ (with $I_{02}$ and $C_2$) and geometrical inductances $L_1$ and $L_2$. The device is controlled with a current bias $I_b$ and a flux current $I_f$ which generates flux $\Phi_a$ through mutual inductance $M$. Transitions can be induced by a microwave current [$I_{r\!f}$]{}, which is coupled to $J1$ via [$C_{r\!f}$]{}. (b) When biased appropriately, the dynamics of the phase difference $\gamma_1$ across the qubit junction are analogous to those of a ball in a one-dimensional tilted washboard potential $U$. The metastable state [$\left| n \right>$]{} differs in energy from [$\left| m \right>$]{} by $\hbar {\ensuremath{\omega_{n m}}}$ and tunnels to the voltage state with a rate [$\Gamma_{n}$]{}. (c) The photograph shows a Nb/AlO$_\text{x}$/Nb device. Not seen is an identical SQUID coupled to this device intended for two-qubit experiments; the second SQUID was kept unbiased throughout the course of this work.](Devicev3.eps){width="3.0in"}

Figure \[FDevice\](a) shows the circuit schematic for our dc SQUID phase qubit.[@Palomaki06a] The qubit junction $J1$ (with critical current $I_{01}$ and capacitance $C_1$) is shown on the left. It is isolated from the current bias source $I_b$ by geometrical inductances $L_1$ and $L_2$ and the second junction $J2$ (with $I_{02}$ and $C_2$). In order to independently control the currents in the two arms of the resulting dc SQUID, a current source $I_f$ applies a flux $\Phi_a$ to the SQUID loop through mutual inductance $M$. Good isolation of the qubit junction is obtained when $L_1 / M \gg 1$ and $L_1 / {\ensuremath{\left( L_2 + L_{J2} \right)}} \gg 1$, where $L_{J2}$ is the Josephson inductance of the isolation junction.[@Martinis02a]

For arbitrary values of the bias current $I_b$ and flux current $I_f$, the dynamics of a dc SQUID can be described by 2 degrees of freedom corresponding to the phase differences across each of the junctions. We, however, operate the device by increasing $I_b$ by $\Delta I_b$ while simultaneously increasing $I_f$ by $L_1 \Delta I_b / M$. This nominally keeps the total current through the isolation junction $J2$ near zero, so that the qubit junction current is roughly $I_b$. Furthermore, one obtains a weak dynamical coupling between the junctions by choosing $L_1$ to be large and biasing the SQUID so that the two junctions are well out of resonance with each other.[@Palomaki06a] In this case, the dynamics of the phase difference $\gamma_1$ across the qubit junction are governed to a good approximation[@Mitra06a] by the Hamiltonian of a single current-biased junction,[@Fulton74a; @Leggett87a] $${\ensuremath{\mathcal{H}}}= \frac{4 E_C}{\hbar^2} p_1^2 - E_J {\ensuremath{\left( \cos \gamma_1
         + \frac{I_b - {\ensuremath{I_{r\!f}}}\cos {\ensuremath{\omega_{r\!f}}}t}{I_{01}} \gamma_1 \right)}}.
\label{eqa}$$ Here, $E_C = e^2 / 2 C_1$ and $E_J = I_{01} {\ensuremath{\Phi_0}}/ 2 \pi$ are the charging energy and Josephson coupling energy of the qubit junction, $p_1 = {\ensuremath{\left( {\ensuremath{\Phi_0}}/ 2 \pi \right)}}^2 C_1 \dot{\gamma_1}$ is the momentum conjugate to $\gamma_1$, and [$I_{r\!f}$]{} is the amplitude of a microwave drive current of frequency [$\omega_{r\!f}$]{}. In quantizing ${\ensuremath{\mathcal{H}}}$, $\gamma_1$ and $p_1$ become operators with ${\ensuremath{\left[ \gamma_1, p_1 \right]}} = i \hbar$.

The second term on the right-hand side of [Eq. [(\[eqa\])]{}]{} defines a one-dimensional tilted washboard potential $U$, sketched in Fig. \[FDevice\](b). Each well is characterized by the classical plasma frequency $\omega_p$ and barrier height $\Delta U$, which can also be expressed as the dimensionless quantity $N_s = \Delta U / \hbar \omega_p$. A single potential well supports roughly $N_s$ metastable states [$\left| n \right>$]{}, where the ground state [$\left| 0 \right>$]{} and first excited state [$\left| 1 \right>$]{} serve as the basis for quantum computation.[@Ramos01a] Motivated by this, the Hamiltonian can be expressed in a discrete representation as $${\ensuremath{\mathcal{H}}}_N = \sum_{n=0}^{N-1} \hbar {\ensuremath{\omega_{0 n}}} {\ensuremath{\left| n \right>}} {\ensuremath{\left< n \right|}}
         + \sum_{n, m = 0}^{N-1} \hbar {\ensuremath{\Omega_{n m}}} {\ensuremath{\left| n \right>}} {\ensuremath{\left< m \right|}} \cos {\ensuremath{\omega_{r\!f}}}t,
\label{eqc}$$ where $N$ is the number of states in the well being considered, $\hbar {\ensuremath{\omega_{n m}}}$ is the energy-level spacing between states [$\left| n \right>$]{} and [$\left| m \right>$]{} (tunable through $I_b$), and $${\ensuremath{\Omega_{n m}}} = \frac{{\ensuremath{\Phi_0}}}{2 \pi} \frac{{\ensuremath{I_{r\!f}}}}{\hbar} {\ensuremath{\left< n \left| \, \gamma_1 \, \right|
                                    m \right>}}
\label{eqf}$$ is a bare Rabi frequency.

The simplified Hamiltonian in [Eq. [(\[eqc\])]{}]{} neglects tunneling through the potential barrier,[@Larkin86a; @Leggett87a; @Chow88a; @Kopietz88a] a process corresponding to the SQUID spontaneously switching from the supercurrent to the voltage state; tunneling from [$\left| n \right>$]{} occurs at an escape rate [$\Gamma_{n}$]{}. As ${\ensuremath{\Gamma_{n+1}}} / {\ensuremath{\Gamma_{n}}} \sim 10^3$ for values of $I_b$ that we studied, tunneling can be exploited to perform state readout. Note also typically ${\ensuremath{\Gamma_{n}}} \ll {\ensuremath{\omega_{n, n+1}}}$ for the lowest few levels; for states near the barrier, this is not true. The total escape rate of the system to the voltage state is $$\Gamma = \frac{1}{\rho_{tot}} \sum_{n=0}^{N-1} \rho_{nn} {\ensuremath{\Gamma_{n}}},
\label{eqg}$$ where $\rho_{nn}$ is the occupation probability of state [$\left| n \right>$]{} and $\rho_{tot} = \sum \rho_{nn}$ is the probability of the system being in the supercurrent state; all of the quantities in the equation are time dependent during a bias ramp.

This description of the qubit makes several simplifications. For one, we have ignored the quantum states of the isolation junction $J2$.[@Mitra06a] For typical bias conditions, we keep the current through this junction small and its first excited state is two to three times higher in energy than that of the qubit. However, the higher excited states of the qubit junction may come into resonance with states of the isolation junction. This approach further ignores the role of $J2$ in determining the bias conditions. Quantum mechanical simulations of the SQUID show that the current through the qubit junction weakly depends on its quantum state;[@Mitra06a] shifts in the qubit current are predicted to be less than 5 nA (much smaller than $I_{01} \approx 18\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$) for our situation, so we have not considered this correction below. As a result, the junction parameters in [Eq. [(\[eqa\])]{}]{} need not be equal to the actual values for the qubit junction.

The simulations described in the sections that follow require values of [$\Gamma_{n}$]{}, [$\omega_{n m}$]{}, and [$\left< n \left| \, \gamma_1 \, \right|
                                    m \right>$]{}, all of which are specified by $I_b$, $I_{01}$, and $C_1$. In the absence of dissipation, analytical expressions for these quantities can be obtained by applying perturbation theory to the cubic approximation of the tilted washboard potential.[@Strauch04a] For the escape rates and energy levels, we instead use numerical solutions to the exact potential, which are more reliable for states near the top of the barrier. These solutions come from two methods that give consistent results: (i) solving the Schrödinger equation with transmission boundary conditions[@Shibata91a] and (ii) complex scaling.[@Moiseyev98a]

We performed experiments on a dc SQUID phase qubit fabricated by Hypres, Inc.[@NIST] using a $30\ {\ensuremath{{\ensuremath{\mathrm{A}}} / {\ensuremath{\mathrm{cm}}}^2}}$ Nb/AlO$_\text{x}$/Nb trilayer process \[see Fig. \[FDevice\](c)\]. The qubit and isolation junctions had areas of $10 \times 10~{\ensuremath{\mu {\ensuremath{\mathrm{m}}}}}^2$ and $7 \times 7~{\ensuremath{\mu {\ensuremath{\mathrm{m}}}}}^2$, respectively, and the inductances were roughly $L_1 = 3.4\ {\ensuremath{{\ensuremath{\mathrm{nH}}}}}$, $L_2 = 30\ {\ensuremath{{\ensuremath{\mathrm{pH}}}}}$, and $M = 13.4\ {\ensuremath{{\ensuremath{\mathrm{pH}}}}}$.[@Palomaki07a] The error in these and other fit parameters reported in this paper is about one unit in the least significant digit. We have not attempted a thorough analysis of the uncertainties in all the parameters due to the complexities of the nonlinear functions involved. The microwave current [$I_{r\!f}$]{} was carried inside the refrigerator on a single length of lossy stainless coax. It was coupled to the qubit via an on-chip $1\ {\ensuremath{{\ensuremath{\mathrm{fF}}}}}$ capacitor [$C_{r\!f}$]{} \[see Fig. \[FDevice\](a)\]; while the large impedance mismatch this produced potentially improved the isolation, it did make an independent calibration of the power reaching the qubit difficult. The device was mounted in a superconducting aluminum box that was attached to the mixing chamber of a dilution refrigerator with a base temperature of 20 mK. The refrigerator, located in an rf shielded room, was surrounded by a mu-metal cylinder. In addition, the measurement and bias lines were protected from noise at room temperature by discrete $LC$ filters and copper powder filters at the mixing chamber.

The escape rate of the qubit junction was measured by simultaneously increasing $I_b$ and $I_f$ as described above, while monitoring the voltage on the current bias line of the SQUID. We measured the time interval between the start of the ramps and when the system tunneled to the voltage state. Repeating this procedure many times ($\sim 10^5$) at a rate of about 250 Hz yielded a histogram of switching times, from which the escape rate could be calculated.[@Fulton74a] Decreasing the repetition rate did not yield a significant change in the escape rate, suggesting that heating due to the device switching to the voltage state had a minimal impact on the measurements. Because the large loop inductance of the SQUID resulted in about 20 possible flux states, the device was initialized to the zero trapped flux state with a flux shaking procedure before each repetition of the measurement cycle.[@Palomaki06a] We typically applied 50 oscillations of the flux current, which yielded a success rate better than 98%.

![\[FEnergyGamma\] Qubit energy-level spectroscopy and tunneling escape rates. (a) Open circles show the resonance frequency of the transition between the ground and first excited states of the qubit, measured at 20 mK. The scatter in the values is indicative of the uncertainty in the measurement. Also plotted are theoretical values of [$\omega_{0 1}$]{} (solid), [$\omega_{1 2}$]{} (dashed), and [$\omega_{2 3}$]{} (dotted) for $I_{01} = 17.930\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ and $C_1 = 4.50\ {\ensuremath{{\ensuremath{\mathrm{pF}}}}}$. (b) At high bias, the measured background escape rate (open circles) agrees with the predicted ground-state escape rate [$\Gamma_{0}$]{} (solid line) for the junction parameters given above. Calculated [$\Gamma_{1}$]{}, [$\Gamma_{2}$]{}, and [$\Gamma_{3}$]{} are plotted as dashed, dotted, and dashed-dotted lines. In both plots, the bottom axes show the total bias current $I_b$, while the top axes indicate the normalized barrier height $N_s$, calculated using the extracted junction parameters.](EnergyGammav3.eps){width="3.0in"}

Figure \[FEnergyGamma\](a) shows a spectrum of transitions for the qubit junction. With a microwave drive of fixed frequency applied to the device, the total escape rate was measured while ramping the biases. The open circles indicate the values of $I_b$ where the microwave drive caused the largest enhancement in the escape rate (over its background value in the absence of microwaves).

As shown by the solid line, we fit these points to the theoretical values of [$\omega_{0 1}$]{}, yielding $I_{01} = 17.930\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ and $C_1 = 4.50\ {\ensuremath{{\ensuremath{\mathrm{pF}}}}}$. These values are likely to differ from the actual critical current and capacitance of the qubit junction, both in principle, due to the simplifications in the Hamiltonian mentioned above, and in practice, as this fitting procedure is sensitive to inaccuracies in the simultaneous ramping of $I_b$ and $I_f$, constant offset flux that biases the SQUID (from, for example, trapped flux near the device), and drifts in the detection electronics. Nonetheless, when we repeated the measurement at elevated temperature, the location of the transitions between higher levels agreed to within 0.2% of predictions for [$\omega_{1 2}$]{} (dashed) and [$\omega_{2 3}$]{} (dotted) obtained with the effective parameters.[@Strauch06a] In addition, two-photon transitions between these levels occurred at high power at the expected frequencies. Multiphoton transitions between the ground and first excited states of a variety of junction qubits have previously been observed and characterized.[@Nakamura01a; @Wallraff03a; @Saito04a] Over the course of three months (during which the refrigerator remained near its base temperature), the junction parameters found from the spectral fits varied by roughly 1%. Thus, we repeated the spectroscopy measurement often to obtain the values of $I_{01}$ and $C_1$ needed for the simulations; these values are listed in the caption of each figure.

The open circles of Fig. \[FEnergyGamma\](b) show the measured background escape rate in the absence of microwaves. The solid line indicates the theoretical value of [$\Gamma_{0}$]{}, calculated with the values of $I_{01}$ and $C_1$ extracted from the spectrum. At low bias current, the measured escape rate exceeds [$\Gamma_{0}$]{}, suggesting the presence of excited state population even though the refrigerator was at 20 mK. Based on experiments described elsewhere,[@Palomaki07a] we believe that features such as the peak near $I_b = 17.73\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ are due to population in [$\left| 2 \right>$]{} generated when noise on the leads at the resonant frequency of the isolation junction matches the [$0 \rightarrow 2$]{} transition frequency of the qubit. At high bias, the escape rates exceed the excitation rates of the noise and the total escape rate converges to $\Gamma_0$; a similar effect occurs for thermal excitations.[@Dutta04a] While noise complicates the situation, the overall behavior suggests that both [$\omega_{0 1}$]{} and [$\Gamma_{0}$]{} of the dc SQUID are described by the same one-dimensional tilted washboard potential. Also shown in Fig. \[FEnergyGamma\](b) are numerical predictions for [$\Gamma_{1}$]{} (dashed), [$\Gamma_{2}$]{} (dotted), and [$\Gamma_{3}$]{} (dashed-dotted).

We observed Rabi oscillations between the ground and first excited states by turning the microwave current [$I_{r\!f}$]{} on when $I_b$ was at the value where ${\ensuremath{\omega_{0 1}}} = {\ensuremath{\omega_{r\!f}}}$ and measuring the time-resolved escape rate while [$I_{r\!f}$]{} remained at a constant value. This serves as a simple method of monitoring the evolution of the state populations.[@Yu02a] Although the biases continued to increase during the Rabi oscillations, the ramp rates were reduced (with $d I_b / dt \approx 0.01\ {\ensuremath{{\ensuremath{\mathrm{A}}} / {\ensuremath{\mathrm{s}}}}}$) before [$I_{r\!f}$]{} was turned on, so that the level spacing ${\ensuremath{\omega_{0 1}}}$ changed by a negligible amount during the escape rate measurement. The symbols in Fig. \[F031806DGtot\] show $\Gamma$ due to a 6.2 GHz microwave drive for a range of powers ${\ensuremath{P_S}}$ at the microwave generator. As expected, the oscillation frequency increases with power and decoherence causes the amplitude of the oscillations to decay with time.

However, there are two unexpected features in the data. First, the escape rate [$\Gamma_{\infty}$]{} at long time (once the oscillations have decayed away) increases with [$P_S$]{}, whereas for an ideal two-level system, the excited state population saturates at high power. Furthermore, the highest measured escape rates far exceed the value of $2.2 \times 10^6\ {\ensuremath{{\ensuremath{\mathrm{1/s}}}}}$ predicted for [$\Gamma_{1}$]{} at this bias current, strongly suggesting that the states [$\left| 2 \right>$]{} and higher are becoming occupied. As the escape rates from these levels are very large, only a small population would be required to account for the observed $\Gamma$. Second, at high power, the oscillation maxima increase over the first few cycles \[see Fig. \[F031806DGtot\](a)\]; this is due to the rise time of the microwave current amplitude, which will be discussed in Sec. \[SDecoherence\].

In Sec. \[SStrong\], we will show how more complete measurements reveal that the higher excited states impact the Rabi oscillation frequencies in ways that are in quantitative agreement with predictions from the Hamiltonian of [Eq. [(\[eqc\])]{}]{}, but are difficult to see in Fig. \[F031806DGtot\]. As shown by the solid lines in the figure, most of the features of the measured escape rate are captured by a multilevel density-matrix simulation which will be described in Sec. \[SDecoherence\].

![\[F031806DGtot\] Rabi oscillations in the escape rate $\Gamma$ were induced at $I_b = 17.746\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ by switching on a microwave current at $t = 0$ with a frequency of 6.2 GHz (resonant with the [$0 \rightarrow 1$]{} transition) and source powers [$P_S$]{} between $-12$ and $-32$ dBm, as labeled. The measurements were taken at 20 mK. The solid lines are from a five-level density-matrix simulation with $I_{01} = 17.930\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$, $C_1 = 4.50\ {\ensuremath{{\ensuremath{\mathrm{pF}}}}}$, $T_1 = 17\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$, and $T_\phi = 16\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$.](031806DGtotv3.eps){width="3.0in"}

Detuning and Strong Field Effects {#SStrong}
=================================

We first model Rabi oscillations using the rotating wave approximation, in the absence of tunneling and dissipation. This provides a way of predicting the oscillation frequency for a wide range of experimental parameters.[@Steffen03a; @Claudon04a; @Meier05a; @Amin06a; @Strauch06a; @Claudon07b] In the rotating frame corresponding to the drive frequency ${\ensuremath{\omega_{r\!f}}}\approx {\ensuremath{\omega_{0 1}}}$, the Hamiltonian of [Eq. [(\[eqc\])]{}]{} simplifies to[@Strauch06a] $${\ensuremath{\mathcal{H}}}_N^{RW} = \sum_{n=0}^{N-1} \hbar \Delta_n {\ensuremath{\left| n \right>}} {\ensuremath{\left< n \right|}}
             + \frac{1}{2} \sum_{\substack{n,m=0 \\ n \neq m}}^{N-1}
               \hbar \Omega_{nm}^\prime {\ensuremath{\left| n \right>}} {\ensuremath{\left< m \right|}},
\label{eqd}$$ where $\Delta_n = {\ensuremath{\omega_{0 n}}} - n {\ensuremath{\omega_{r\!f}}}$ and $$\Omega_{nm}^\prime = {\ensuremath{\Omega_{n m}}} \sum_{s=0}^1 J_{m-n+2s-1}
                       {\ensuremath{\left( \frac{{\ensuremath{\Omega_{n n}}} - {\ensuremath{\Omega_{m m}}}}{{\ensuremath{\omega_{r\!f}}}} \right)}}$$ for $n < m$ and $\Omega_{mn}^\prime = \Omega_{nm}^\prime$. Here, $J_n {\ensuremath{\left( x \right)}}$ is the $n$th order Bessel function and $\Omega_{nm}$ is defined by [Eq. [(\[eqf\])]{}]{}.

Differences between the $N$ eigenvalues of this Hamiltonian specify effective multilevel Rabi frequencies or modes of the system. We find it convenient to label these frequency differences by the states [$\left| n \right>$]{} with the largest weight in the two corresponding eigenfunctions. For the parameter regime of interest here, the differences can be uniquely classified by two states, which we denote by [$\Omega_{R,n m}$]{}. For example, [$\Omega_{R,0 2}$]{} denotes the Rabi oscillation frequency between eigenstates that describe a two-photon transition between [$\left| 0 \right>$]{} and [$\left| 2 \right>$]{}. While approximate analytical solutions can be obtained for a three-level junction system,[@Strauch06a] we numerically found the eigenvalues of the rotating wave Hamiltonian for systems with up to seven levels in the simulations that follow.

![\[F031806Dfosc\] Rabi oscillation frequency [$\Omega_{R,0 1}$]{} at fixed bias as a function of microwave current [$I_{r\!f}$]{}. Extracted values from data (including the plots in Fig. \[F031806DGtot\]) are shown as circles, while the rotating wave solution is shown for two- (dashed line) and five- (solid) level simulations, calculated using $I_{01} = 17.930\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ and $C_1 = 4.50\ {\ensuremath{{\ensuremath{\mathrm{pF}}}}}$ with ${\ensuremath{\omega_{r\!f}}}/ 2 \pi = 6.2\ {\ensuremath{{\ensuremath{\mathrm{GHz}}}}}$.](031806Dfoscv3.eps){width="3.0in"}

We fit the escape rates in Fig. \[F031806DGtot\] (and additional data for other powers not shown) to a decaying sinusoid with an offset. The extracted frequencies are shown with circles in Fig. \[F031806Dfosc\]. To compare to theory, [$\Omega_{R,0 1}$]{}, calculated using the rotating wave solution for a system with five levels, is shown with a solid line. The implied assumption that the oscillation frequencies of $\Gamma$ and $\rho_{11}$ are equal, even at high power in a multilevel system, will be addressed in Sec. \[SSummary\]. In plotting the data, we have introduced a single fitting parameter $117\ {\ensuremath{{\ensuremath{{\ensuremath{\mathrm{nA}}}}}/ \sqrt{{\ensuremath{{\ensuremath{\mathrm{mW}}}}}}}}$ that converts the power [$P_S$]{} at the microwave source to the current amplitude [$I_{r\!f}$]{} at the qubit. Good agreement is found over the full range of power.

As [$I_{r\!f}$]{} increases in Fig. \[F031806Dfosc\], the oscillation frequency is smaller than the expected linear relationship for a two-level system (dashed line). This effect is a hallmark of a multilevel system and has been previously observed in a similar phase qubit.[@Claudon04a; @Claudon07b] There are two distinct phenomena that affect [$0 \rightarrow 1$]{} Rabi oscillations in such a device.[@Goorden03a; @Meier05a; @Strauch06a; @Shevchenko07a] To describe these, we must first define resonance as occurring when the Rabi frequency [$\Omega_{R,0 1}$]{} is at a minimum, as the detuning between the microwave drive frequency [$\omega_{r\!f}$]{} and level spacing [$\omega_{0 1}$]{} is varied (for a fixed drive power). In a two-level system, this happens for ${\ensuremath{\omega_{0 1}}} = {\ensuremath{\omega_{r\!f}}}$. For our phase qubit, resonance occurs when ${\ensuremath{\omega_{0 1}}} < {\ensuremath{\omega_{r\!f}}}$ due to the decreasing level spacings with increasing state [$\left| n \right>$]{} \[see Fig. \[FEnergyGamma\](a)\]. This shift in resonance is an analog of the ac Stark effect. Resonance shifts occur under strong driving in other superconducting qubits as well.[@Schuster05a; @Dutton06a]

However, this shift does not explain the data shown in Fig. \[F031806Dfosc\], which were taken at a fixed bias current (and thus off resonance at high power). The higher levels also affect the frequency of the Rabi oscillations. We will refer to the minimum value of [$\Omega_{R,n m}$]{} as the on-resonance Rabi frequency [$\Omega^{min}_{R,n m}$]{}. The suppression of [$\Omega^{min}_{R,0 1}$]{} below [$\Omega_{0 1}$]{} leads to the effect seen in Fig. \[F031806Dfosc\].

Both of these effects become significant as [$\Omega_{0 1}$]{} approaches ${\ensuremath{\omega_{0 1}}} - {\ensuremath{\omega_{1 2}}}$, which is a measure of the anharmonicity of the system;[@Meier05a; @Amin06a; @Strauch06a] in the case of Fig. \[F031806Dfosc\], ${\ensuremath{\omega_{0 1}}} / 2 \pi = 6.2\ {\ensuremath{{\ensuremath{\mathrm{GHz}}}}}$ and ${\ensuremath{\omega_{1 2}}} / 2 \pi = 5.5\ {\ensuremath{{\ensuremath{\mathrm{GHz}}}}}$. Clearly, these shifts need to be considered when working at high power or at low current bias.

In order to follow experimentally the shift of the resonance condition, it was necessary to measure Rabi oscillations for different detunings of the microwave drive. We chose to do this by keeping the drive frequency [$\omega_{r\!f}$]{} fixed and changing the level spacing [$\omega_{0 1}$]{} (through $I_b$), because the power transmitted by the microwave lines had a nontrivial frequency dependence. Figure \[F010206H1\](a) shows a grayscale plot of Rabi oscillations measured from such an experiment, where black represents a high escape rate. Each horizontal line is the escape rate versus time due to a microwave current of 6.5 GHz and $-11$ dBm, which was turned on at the value of the current bias $I_b$ indicated on the vertical axis. While the measurements were performed at 110 mK, this is not expected to have a significant impact on the Rabi oscillations, as the temperature was well below $\hbar {\ensuremath{\omega_{0 1}}} / k_B \approx 325\ {\ensuremath{{\ensuremath{\mathrm{mK}}}}}$.[@Lisenfeld07a]

![\[F010206H1\] (Color online) Multiphoton, multilevel Rabi oscillations plotted in the time and frequency domains. (a) The escape rate $\Gamma$ (measured at 110 mK) is plotted as a function of the time after which a 6.5 GHz, $-11$ dBm microwave drive was turned on and the current bias $I_b$ of the qubit; $\Gamma$ ranges from 0 (white) to $3 \times 10^8\ {\ensuremath{{\ensuremath{\mathrm{1/s}}}}}$ (black). (b) The normalized power spectral density of the time-domain data from $t = 1$ to 45 ns is shown with a grayscale plot. The dashed line segments indicate the Rabi frequencies obtained from the rotating wave model for transitions involving (from top to bottom) 1, 2, 3, and 4 photons, evaluated with junction parameters $I_{01} = 17.828\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ and $C_1 = 4.52\ {\ensuremath{{\ensuremath{\mathrm{pF}}}}}$, and microwave current ${\ensuremath{I_{r\!f}}}= 24.4\ {\ensuremath{{\ensuremath{\mathrm{nA}}}}}$. Corresponding grayscale plots calculated with a seven-level density-matrix simulation are shown in (c) and (d).](010206H1v3.eps){width="3.0in"}

From Fig. \[F010206H1\](a), we see that the oscillation frequency depends on the current bias. This variation can be seen more readily in Fig. \[F010206H1\](b), which shows the power spectral density (calculated as the absolute square of the discrete Fourier transform) of the escape rate data in Fig. \[F010206H1\](a) between $t = 1$ and $45\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$. For this plot, each horizontal line has been normalized to its maximum value (black) in order to emphasize the location of the dominant frequency. Three distinct bands are visible.

For this data set, the level spacing ${\ensuremath{\omega_{0 1}}} / 2 \pi$ is equal to the microwave frequency ${\ensuremath{\omega_{r\!f}}}/ 2 \pi = 6.5\ {\ensuremath{{\ensuremath{\mathrm{GHz}}}}}$ at $I_b = 17.614\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$. The band with the highest current in Fig. \[F010206H1\](b) is centered about $I_b = 17.624\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$, suggesting that [$0 \rightarrow 1$]{} Rabi oscillations are the dominant process near this bias. For slightly higher or lower $I_b$, the oscillation frequency increases as ${\ensuremath{\Omega_{R,0 1}}} \approx \sqrt{\Omega_{01}^{\prime2} + {\ensuremath{\left( {\ensuremath{\omega_{r\!f}}}- {\ensuremath{\omega_{0 1}}} \right)}}^2}$, in agreement with simple two-level Rabi theory, leading to the curved band in the grayscale plot.

The other bands correspond to oscillations between the ground state and higher excited states. Due to the anharmonic level structure, ${\ensuremath{\omega_{0 2}}} / 4 \pi$ is 6.5 GHz at a smaller current bias $I_b = 17.594\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$; thus a second band appears there, corresponding to two-photon [$0 \rightarrow 2$]{} Rabi oscillations. Similarly, three-photon [$0 \rightarrow 3$]{} oscillations are visible near $17.572\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ where ${\ensuremath{\omega_{0 3}}} / 6 \pi = 6.5\ {\ensuremath{{\ensuremath{\mathrm{GHz}}}}}$. Finally, large escape rates occur near $17.549\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ corresponding to a four-photon [$0 \rightarrow 4$]{} transition, although no oscillations are apparent in Fig. \[F010206H1\](a). Note that the effective junction parameters $I_{01}$ and $C_1$ (given in the caption) used to predict [$\omega_{n m}$]{} and other level properties are slightly different than those for Figs. \[FEnergyGamma\]–\[F031806Dfosc\], as the data sets were taken two months apart.

The rotating wave solution provides a simple way to predict the oscillation frequencies. Calculations of [$\Omega_{R,0 n}$]{} using [Eq. [(\[eqd\])]{}]{} for a seven-level system are shown as dashed curves in Fig. \[F010206H1\](b) for the four lowest multiphoton transitions ($n = 1, 2, 3, 4$). The microwave amplitude ${\ensuremath{I_{r\!f}}}= 24.4\ {\ensuremath{{\ensuremath{\mathrm{nA}}}}}$ is the only free parameter in the calculation and the rotating wave solution using this value reproduces the oscillation frequencies of the different processes well, even at large detuning.

Figure \[F010206H1\](b) shows that the minimum oscillation frequency ${\ensuremath{\Omega^{min}_{R,0 1}}} / 2 \pi = 540\ {\ensuremath{{\ensuremath{\mathrm{MHz}}}}}$ of the first (experimental) band occurs at $I_b = 17.624\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$, for which ${\ensuremath{\omega_{0 1}}} / 2 \pi = 6.4\ {\ensuremath{{\ensuremath{\mathrm{GHz}}}}}$. This again indicates an ac Stark shift of this transition, which we denote by $\Delta {\ensuremath{\omega_{0 1}}} \equiv {\ensuremath{\omega_{r\!f}}}- {\ensuremath{\omega_{0 1}}} \approx 2 \pi \times 100\ {\ensuremath{{\ensuremath{\mathrm{MHz}}}}}$. In addition, the higher levels have suppressed the oscillation frequency below the bare Rabi frequency of ${\ensuremath{\Omega_{0 1}}} / 2 \pi = 620\ {\ensuremath{{\ensuremath{\mathrm{MHz}}}}}$ \[calculated with [Eq. [(\[eqf\])]{}]{}\].

![\[F032206MNstats\] (Color online) The (a) on-resonance Rabi oscillation frequencies [$\Omega^{min}_{R,0 1}$]{} and [$\Omega^{min}_{R,0 2}$]{} and (b) resonance frequency shifts $\Delta {\ensuremath{\omega_{0 1}}} = {\ensuremath{\omega_{r\!f}}}- {\ensuremath{\omega_{0 1}}}$ and $\Delta {\ensuremath{\omega_{0 2}}} = 2 {\ensuremath{\omega_{r\!f}}}- {\ensuremath{\omega_{0 2}}}$ are plotted as a function of the microwave current, for data taken at $110\ {\ensuremath{{\ensuremath{\mathrm{mK}}}}}$ with a microwave drive of frequency ${\ensuremath{\omega_{r\!f}}}/ 2 \pi = 6.5\ {\ensuremath{{\ensuremath{\mathrm{GHz}}}}}$ and powers ${\ensuremath{P_S}}= -23, -20, -17, -15, -10\ {\ensuremath{{\ensuremath{\mathrm{dBm}}}}}$. Values extracted from data for the [$0 \rightarrow 1$]{} ([$0 \rightarrow 2$]{}) transition are plotted as open circles (filled squares), while five-level rotating wave solutions for a junction with $I_{01} = 17.736\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ and $C_1 = 4.49\ {\ensuremath{{\ensuremath{\mathrm{pF}}}}}$ are shown as solid (dashed) lines. In (a), the dotted line is from a simulation of a two-level system.](032206MNstatsv3.eps){width="3.0in"}

We repeated this analysis of resonance for five different microwave powers, with data taken at a later date. Figure \[F032206MNstats\] shows experimental results for the on-resonance Rabi frequencies ${\ensuremath{\Omega^{min}_{R,n m}}}$ and Stark shifts $\Delta {\ensuremath{\omega_{n m}}}$ with the corresponding results from a five-level rotating wave solution for the [$0 \rightarrow 1$]{} (circles for data and solid lines for theory) and two-photon [$0 \rightarrow 2$]{} (squares and dashed lines) transitions. Here, the power calibration is $84\ {\ensuremath{{\ensuremath{{\ensuremath{\mathrm{nA}}}}}/ \sqrt{{\ensuremath{{\ensuremath{\mathrm{mW}}}}}}}}$ for ${\ensuremath{\omega_{r\!f}}}/ 2 \pi = 6.5\ {\ensuremath{{\ensuremath{\mathrm{GHz}}}}}$. Figure \[F032206MNstats\](a) differs from Fig. \[F031806Dfosc\], because in the former, $I_b$ was varied at each power to give the minimum oscillation frequency; by staying on resonance in this way, the effect of the higher levels on the [$0 \rightarrow 1$]{} oscillations is maximized.

The resonant oscillation frequencies [$\Omega^{min}_{R,0 1}$]{} and [$\Omega^{min}_{R,0 2}$]{} in Fig. \[F032206MNstats\](a) are well described by [Eq. [(\[eqd\])]{}]{} over the full range of [$I_{r\!f}$]{}. The deviation between the [$0 \rightarrow 1$]{} oscillation frequency and the values expected in a two-level system (dotted line) increases with [$I_{r\!f}$]{}. Similar measurements taken at a higher $I_b$ show a smaller frequency suppression over a similar range of [$I_{r\!f}$]{},[@Strauch06a] as expected for a system with stronger anharmonicity.

The resonance shifts $\Delta {\ensuremath{\omega_{0 1}}}$ and $\Delta {\ensuremath{\omega_{0 2}}}$ in Fig. \[F032206MNstats\](b), however, differ significantly from the model predictions. As a change of 10 MHz in ${\ensuremath{\omega_{0 1}}} / 2 \pi$ corresponds to a roughly 1 nA change in $I_b$, the discrepancy is difficult to see in Fig. \[F010206H1\](b). The indicated uncertainty in the experimental points in Fig. \[F032206MNstats\](b) is roughly 5 MHz due to errors in the calibration of ${\ensuremath{\omega_{0 1}}} {\ensuremath{\left( I_b \right)}}$. In addition the uncertainty is somewhat larger at low power, where the relatively small total escape rates result in poor counting statistics, and high power, where the weak dependence on detuning makes it difficult to identify the resonant level spacing. While the general trend of the shifts is consistent with the model, the scatter in the data is large and further work would be needed to determine if there are true deviations from the multilevel theory.

Decoherence {#SDecoherence}
===========

We now model the time dependence of the escape rate measured in the experiment. In order to do this, several additions have to be made to the treatment of the system given in Sec. \[SStrong\]: the effects of tunneling, other sources of damping and noise, and experimental limitations. The density-matrix formalism[@Smith78a; @Blum96a] provides a straightforward scheme for including nonunitary processes and has been previously applied to the lowest two or three levels of the phase qubit.[@Yu02a; @Kosugi05a; @Meier05a; @Amin06a; @Claudon07b; @Shevchenko07a] For a system with $N$ levels, we assume that the evolution of the qubit’s reduced density matrix [$\rho$]{} is given by the modified Liouville–von Neumann equation $${\ensuremath{\frac{\partial {\ensuremath{\rho}}}{\partial t}}} = -\frac{i}{\hbar} {\ensuremath{\left[ {\ensuremath{\mathcal{H}}}_N, {\ensuremath{\rho}}\right]}} - G {\ensuremath{\rho}}- R {\ensuremath{\rho}}- D {\ensuremath{\rho}}.
\label{eqb}$$ Here, we use the discrete Hamiltonian in [Eq. [(\[eqc\])]{}]{}, so that the diagonal elements $\rho_{nn}$ give the occupancy of the states [$\left| n \right>$]{}. Tunneling is characterized by the tensor $G$, where ${\ensuremath{\left[ G {\ensuremath{\rho}}\right]}}_{nm} = {\ensuremath{\left( {\ensuremath{\Gamma_{n}}} + {\ensuremath{\Gamma_{m}}} \right)}} \rho_{nm} / 2$, which leads to a decay of all elements of [$\rho$]{}. The tensors $R$ and $D$ account for two distinct decoherence mechanisms.

The Bloch-Redfield tensor $R$ models the effects of the system being in equilibrium with a thermal bath at temperature $T$.[@Burkard04a] This tensor leads to the decay of the diagonal elements of $\rho$ (dissipation) as well as the off-diagonal elements (decoherence). The coupling to the bath, assumed linear in $\gamma_1$, is parametrized by $R_1 {\ensuremath{\left( \omega \right)}}$, which is the inverse of the real part of the total admittance that shunts the qubit junction, evaluated at angular frequency $\omega$. In [Eq. [(\[eqb\])]{}]{}, $${\ensuremath{\left[ R {\ensuremath{\rho}}\right]}}_{nm} = \sum_{k,l=0}^{N-1} R_{nmkl} \rho_{kl},
\label{eqe}$$ where $$R_{nmkl} = - \gamma_{lmnk} - \gamma_{knml}
             + \delta_{lm} \sum_{r=0}^{N-1} \gamma_{nrrk}
             + \delta_{nk} \sum_{r=0}^{N-1} \gamma_{mrrl}$$ and

$$\gamma_{lmnk} = \frac{1}{2 \hbar} {\ensuremath{\left( \frac{{\ensuremath{\Phi_0}}}{2 \pi} \right)}}^2
      \frac{{\ensuremath{\left< l \left| \, \gamma_1 \, \right|
                                    m \right>}} {\ensuremath{\left< n \left| \, \gamma_1 \, \right|
                                    k \right>}}}
      {R_1 {\ensuremath{\left( {\ensuremath{\omega_{n k}}} \right)}}}
      {\ensuremath{\left[ {\ensuremath{\left( 1 - \delta_{nk} \right)}} {\ensuremath{\omega_{n k}}}
      \frac{\exp {\ensuremath{\left( -\frac{\hbar {\ensuremath{\omega_{n k}}} \mathrm{sgn} {\ensuremath{\left( n - k \right)}}}
             {2 k_B T} \right)}}}
           {\sinh {\ensuremath{\left( \hbar {\ensuremath{\omega_{n k}}} / 2 k_B T \right)}}}
      + \delta_{nk} \frac{2 k_B T}{\hbar} \right]}}.$$

In the following, we assume that $R_1$ is independent of frequency, in which case we can define a dissipation time $T_1 = R_1 C_1$. In this limit, [Eq. [(\[eqe\])]{}]{} gives well known interlevel transition rates.[@Larkin86a; @Chow88a; @Xu05b] For example, thermal excitation from [$\left| n \right>$]{} to [$\left| m \right>$]{} (with $m > n$) occurs at a rate $$\begin{aligned}
  W_{mn}^+ & = & 2 \gamma_{nmmn}  \nonumber \\
           & = & \frac{2}{\hbar} {\ensuremath{\left( \frac{{\ensuremath{\Phi_0}}}{2 \pi} \right)}}^2
                 {\ensuremath{\left( \frac{{\ensuremath{\omega_{n m}}}}{R_1} \right)}} \frac{{\ensuremath{\left| {\ensuremath{\left< n \left| \, \gamma_1 \, \right|
                                    m \right>}} \right|}}^2}
                 {\exp {\ensuremath{\left( \hbar {\ensuremath{\omega_{n m}}} / k_B T \right)}} - 1}.\end{aligned}$$ As required by detailed balance, decay from [$\left| m \right>$]{} to [$\left| n \right>$]{} occurs at a rate $W_{nm}^- = W_{mn}^+ \exp {\ensuremath{\left( \hbar {\ensuremath{\omega_{n m}}} / k_B T \right)}}$. In addition, [Eq. [(\[eqe\])]{}]{} specifies the decoherence due to this dissipation. Given ${\ensuremath{\omega_{n m}}}$ and the matrix elements of $\gamma_1$, all of the thermal rates are specified by $T_1$, which is set to $17\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ for the simulations below. This value comes from additional measurements of the escape rate of the device in the absence of microwaves over a range of temperatures.[@Dutta04a]

We find that this treatment of dissipation alone is insufficient to capture the decay of Rabi oscillations, suggesting that an additional decoherence mechanism is present. This we model by the tensor $D$ in [Eq. [(\[eqb\])]{}]{}, which has the form $$D {\ensuremath{\rho}}= \sum_n \lambda_n {\ensuremath{\left( L_n {\ensuremath{\rho}}L_n^\dag
        - \frac{1}{2} L_n^\dag L_n {\ensuremath{\rho}}- \frac{1}{2} {\ensuremath{\rho}}L_n^\dag L_n \right)}},
\label{eqh}$$ where $L_n$ are Lindblad operators with strengths $\lambda_n$. The best overall agreement with the measurements (see further discussion below) is found for a set of operators $L_n = {\ensuremath{\left| n \right>}} {\ensuremath{\left< n \right|}}$, where $\lambda_n = 1 / T_\phi$ and $n$ ranges from 0 to $N-1$. This leads to an exponential decay of each of the off-diagonal elements of [$\rho$]{} with a common time constant given by the dephasing time $T_\phi$ and no change in the diagonal elements. In the simulations of this section, we set $T_\phi = 16\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$, yielding a coherence time $T_2 = {\ensuremath{\left[ 1 / {\ensuremath{\left( 2 T_1 \right)}} + 1 / T_\phi \right]}}^{-1} = 10.9\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$.

For $N = 2$, [Eq. [(\[eqb\])]{}]{} reduces to the optical Bloch equations,[@Smith78a; @Blum96a] for which analytical solutions exist. However, for $N > 2$ the master equation is easily numerically integrated to obtain the time dependence of [$\rho$]{}, without making the rotating wave approximation. With the state occupation probabilities $\rho_{nn}$ in hand, the total escape rate can be calculated with [Eq. [(\[eqg\])]{}]{} and compared with experiment.

As mentioned earlier, one of the more striking features of the series of [$0 \rightarrow 1$]{} Rabi oscillations in Fig. \[F031806DGtot\] is that [$\Gamma_{\infty}$]{}, which we define as the steady escape rate the system approaches as the oscillations decay away, increases over the full range of measured power. The circles in Fig. \[F031806DGeq\] show experimental values of [$\Gamma_{\infty}$]{} as a function of the bare Rabi frequency [$\Omega_{0 1}$]{}, which was calculated from the microwave source power by using the fit in Fig. \[F031806Dfosc\]. We will now use the density-matrix simulations to understand these escape rates.

![\[F031806DGeq\] Long-time escape rate [$\Gamma_{\infty}$]{} as a function of the bare Rabi frequency [$\Omega_{0 1}$]{}. (a) [$\Gamma_{\infty}$]{} is plotted for a density-matrix simulation of a system with two (dashed-dotted), three (dashed), four (dotted), and five (solid) levels, $I_{01} = 17.930\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$, and $C_1 = 4.50\ {\ensuremath{{\ensuremath{\mathrm{pF}}}}}$. (b) [$\Gamma_{\infty}$]{} at low power (for a five-level system) is affected by the inclusion of a small microwave current at [$\omega_{0 2}$]{} (dashed) and inhomogeneous broadening in addition to this noise current (solid). One set of experimental data, for source powers between $-68$ and $-11\ {\ensuremath{{\ensuremath{\mathrm{dBm}}}}}$ at 6.2 GHz, is plotted in both panels with open circles.](031806DGeqv3.eps){width="3.0in"}

The simulations suggest that [$\Gamma_{\infty}$]{} is not sensitive to the time evolution of [$\rho$]{}, and is only weakly dependent on $T_1$ and $T_\phi$ at high power. However, [$\Gamma_{\infty}$]{} does depend on the number of states that participate in the dynamics and their individual escape rates [$\Gamma_{n}$]{}. It should be noted that the actual populations do not reach steady state, but continuously decay due to tunneling. The dashed-dotted line in Fig. \[F031806DGeq\](a) shows the calculated [$\Gamma_{\infty}$]{} for a two-level system ($N = 2$). It increases up to a value of ${\ensuremath{\Gamma_{1}}} / 2$ at roughly ${\ensuremath{\Omega_{0 1}}} = 1 / \sqrt{T_1 T_2}$; the data show only a subtle shoulder, masked by an overall steady increase, near this escape rate. Results for three (dashed), four (dotted), and five (solid) levels are also plotted in the figure, displaying increasing agreement with experiment.

This agreement between theory and experiment at high power is an indication that transitions to the higher states, or leakage in the context of quantum computation, are occurring as expected. As power increases, the occupation of higher excited states increases and this effect is magnified by their larger escape rates. For example, with the five-level simulation at ${\ensuremath{\Omega_{0 1}}} / 2 \pi = 1\ {\ensuremath{{\ensuremath{\mathrm{GHz}}}}}$ the state [$\left| 3 \right>$]{} has a 2% occupation probability, but accounts for 60% of [$\Gamma_{\infty}$]{}. Note that microwave pulse shaping will not reduce this effect. That is, while leakage can be minimized at the end of a single qubit operation, higher states are always populated during the pulse.[@Steffen03a; @Meier05a; @Amin06a] Such pulses will have enhanced escape rates comparable to those seen here for a constant [$I_{r\!f}$]{}.

These simulations do not do a good job of explaining the value of [$\Gamma_{\infty}$]{} at very low power. The reason for this can be seen in Fig. \[FEnergyGamma\](b), where the measured escape rate in the absence of microwaves exceeds the predicted value of [$\Gamma_{0}$]{}. In particular, at the bias current of $I_b = 17.746\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ where the Rabi oscillations were performed, the measured escape rate is $\sim 2 \times 10^4\ {\ensuremath{{\ensuremath{\mathrm{1/s}}}}}$. While this excess escape rate could be attributed to a thermal population at $T = 56\ {\ensuremath{{\ensuremath{\mathrm{mK}}}}}$ (whereas the refrigerator thermometer read 20 mK), raising the temperature of the simulation leads to a significant enhancement in [$\Gamma_{\infty}$]{} that extends up to moderate microwave power, an effect not seen in the data. Instead, separate measurements[@Palomaki07a] indicate that the spurious features in the measured $\Gamma$ are largely due to a population in [$\left| 2 \right>$]{}, beyond that expected from a thermal bath at 20 mK. Assuming that this was a result of noise on the bias lines, we included in the simulation an additional microwave source at a frequency of [$\omega_{0 2}$]{}; a microwave current amplitude of 0.2 nA reproduces the background escape rate measured at $I_b = 17.746\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ (corresponding to Figs. \[F031806DGtot\], \[F031806Dfosc\], and \[F031806DGeq\]) and leads to a population $\rho_{22} = 3 \times 10^{-5}$. The coherent effects of this drive are negligible, as ${\ensuremath{\Omega_{0 2}}} / 2 \pi = 0.5\ {\ensuremath{{\ensuremath{\mathrm{MHz}}}}}$, which is much smaller than $1/T_1$ or $1/T_\phi$. The results of a five-level simulation with the extra source included are shown with a dashed line in Fig. \[F031806DGeq\](b). Compared to the simpler simulations in Fig. \[F031806DGeq\](a), they show improved agreement at the lowest powers, with little change above ${\ensuremath{\Omega_{0 1}}} / 2 \pi = 10\ {\ensuremath{{\ensuremath{\mathrm{MHz}}}}}$.

To this point, we have not included effects from inhomogeneous broadening, although our spectroscopic measurements suggest its presence. For $T_2 = 10.9\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$, which describes many of the experiments below, the expected full width at half maximum of a resonance peak is roughly 30 MHz. However, we commonly find peak widths of 50 MHz at bias currents where tunneling makes a negligible contribution. The remaining broadening may be due to current noise at frequencies much lower than $1 / T_1$, whose contribution to the width scales with the spectral slope $d {\ensuremath{\omega_{0 1}}} / d I_b$.[@Berkley03b] This can be modeled by taking into account the frequency content of the noise and using the stochastic Bloch equations.[@Xu05b] Instead, we mimicked the inferred spectroscopic broadening simply by running the simulation for a range of bias currents and then convolving the resulting escape rate (at a given time) and a Gaussian with standard deviation $\sigma_I = 1.5\ {\ensuremath{{\ensuremath{\mathrm{nA}}}}}$ (corresponding to a 35 MHz spread in ${\ensuremath{\omega_{0 1}}} / 2 \pi$).

Calculations of [$\Gamma_{\infty}$]{} for a five-level system with noise at [$\omega_{0 2}$]{} and this inhomogeneous broadening are drawn with a solid line in Fig. \[F031806DGeq\](b). The extra broadening has negligible effect at high power (where Rabi oscillations are observed), but does bring the simulation into better agreement with data near ${\ensuremath{\Omega_{0 1}}} / 2 \pi = 10\ {\ensuremath{{\ensuremath{\mathrm{MHz}}}}}$. The small remaining discrepancy could be due to detuning arising from a misidentification of the mean value of $I_b$. For example, performing the simulation at 17.745 rather than $17.746\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ results in an underestimate of the measured values. Nevertheless, the overall agreement over three decades in the Rabi frequency is very good. While the calibration of the microwave current [$I_{r\!f}$]{} came from a fit to data, one conversion factor reproduces both the oscillation frequencies (Fig. \[F031806Dfosc\]) and long-time escape rates (Fig. \[F031806DGeq\]).

![\[F032106FN10\] (Color online) Response to a microwave pulse. (a) The measured escape rate $\Gamma$ (circles), due to a 6.2 GHz, $-20\ {\ensuremath{{\ensuremath{\mathrm{dBm}}}}}$ microwave pulse nominally 30 ns long, shows Rabi oscillations followed by a decay governed by multiple time constants. The solid line is the result of a five-level density-matrix simulation for a junction with $I_{01} = 17.730\ {\ensuremath{\mu {\ensuremath{\mathrm{A}}}}}$ and $C_1 = 4.46\ {\ensuremath{{\ensuremath{\mathrm{pF}}}}}$, using the measured microwave pulse amplitude (with a maximum of 11.7 nA) shown in the inset. The simulation also produces (b) the contribution to the escape rate and (c) the normalized occupation for each of the five levels, as labeled.](032106FN10v3.eps){width="3.0in"}

We next consider the time dependence of the escape rate for the data plotted in Fig. \[F032106FN10\]. Here, a 6.2 GHz microwave pulse nominally 30 ns long was applied on resonance with the [$0 \rightarrow 1$]{} transition of the qubit junction. The measured escape rate shows Rabi oscillations followed by a decay back to the ground state once the microwave drive has turned off. This decay appears to be governed by three time constants. Nontrivial decays have previously been reported in phase qubits[@Martinis02a] and we have found them in several of our devices.

Accurately simulating this experiment requires knowledge of the time dependence of the microwave pulse, which was created by the internal gate of a Hewlett-Packard 83731B synthesized source,[@NIST] without any further filtering. We measured the pulse envelope at the source’s output using a digital sampling oscilloscope \[see inset of Fig. \[F032106FN10\](a)\]. Ignoring any distortion of the pulse before it reached the qubit junction, this was taken as the microwave amplitude ${\ensuremath{I_{r\!f}}}{\ensuremath{\left( t \right)}}$. For the microwave power used in Fig. \[F032106FN10\], the maximum value of [$I_{r\!f}$]{} was taken to be 11.7 nA (a value obtained from the power calibration in Fig. \[F031806Dfosc\]). The solid line in Fig. \[F032106FN10\](a) shows the calculated escape rate for a five-level simulation in the presence of noise at [$\omega_{0 2}$]{} and inhomogeneous broadening. For this time-domain plot and others discussed below, we have convolved the simulation and a Gaussian with a full width of 150 ps to remove very fast, small oscillations (due to highly detuned multiphoton processes) that could not be seen in the experiment due to insufficient time resolution.

The main part of the Rabi oscillation is reproduced well. In particular, the second oscillation maximum has a larger escape rate than the first, due to the 7 ns it takes for [$I_{r\!f}$]{} to reach its maximum. The first part of the decay is also reproduced, which from the simulation should correspond to the emptying of state [$\left| 2 \right>$]{}. However, the data also show a slow decay constant longer than 50 ns that the simulation does not explain. The longer time is inconsistent with our thermal measurement of $T_1$ and instead may be indicative of the qubit interacting with an additional quantum system. It appears, though, that this extra degree of freedom does not significantly affect the description of the Rabi oscillations. Note that if this longer time constant were the dominant relaxation process at high power, two-level saturation would have occurred at a lower power in Fig. \[F031806DGeq\](a).

The density-matrix model can also be used to predict the escape rates for the power series in Fig. \[F031806DGtot\]; calculations are plotted with solid lines in that figure. The presence of the noise signal at ${\ensuremath{\omega_{0 2}}}$ has little effect and inhomogeneous broadening decreases the escape rate at the lowest microwave power slightly. The maximum value of ${\ensuremath{I_{r\!f}}}{\ensuremath{\left( t \right)}}$ was again calculated for each power using the fit in Fig. \[F031806Dfosc\]. The $T_\phi$ value of 16 ns used in the simulations was chosen to best reproduce the decay envelopes over the full range of powers. At the highest power, there is a discrepancy in the oscillation maxima and minima, perhaps indicative of inaccurate parameters for levels [$\left| 3 \right>$]{} and [$\left| 4 \right>$]{}. The envelope of the microwave turn-on is also reflected in the shape of the escape rate.

As a final test of the model, we examine multiphoton transitions in the system. Figure \[F010206H1\](c) shows a grayscale plot of the escape rate calculated with a seven-level density-matrix simulation. Nearly all of the features seen in the data of Fig. \[F010206H1\](a) are present in the simulation. As the gray scales are identical, a small quantitative disagreement is visible, particularly for the [$0 \rightarrow 3$]{} three-photon transition. Figure \[F010206H1\](d) shows the normalized power spectral density calculated from Fig. \[F010206H1\](c); it agrees well with the data in Fig. \[F010206H1\](b).

Figure \[F010206H1lines\] shows line cuts of the time-domain data (circles) and simulation (solid lines) from Fig. \[F010206H1\] at bias currents of (a) 17.623, (b) 17.596, and (c) 17.571 [$\mu {\ensuremath{\mathrm{A}}}$]{}. These values of $I_b$ correspond to the resonances of the one, two, and three-photon transitions at the high power at which the data were taken. While the decay time and the long-time escape rate [$\Gamma_{\infty}$]{} are reproduced well for the three transitions, the first few nanoseconds are not captured fully. This could be due to distortion of the microwave current by the coupling capacitor [$C_{r\!f}$]{} or other line mismatches.

![\[F010206H1lines\] (Color online) Comparison of possible damping scenarios for multi-photon transitions. The measured escape rates from Fig. \[F010206H1\](a) for the (a) [$0 \rightarrow 1$]{}, (b) [$0 \rightarrow 2$]{}, and (c) [$0 \rightarrow 3$]{} transitions are plotted with circles. The solid lines were calculated with the same simulation (with a common dephasing time $T_\phi$) whose results are shown in Fig. \[F010206H1\](c). The dotted lines were calculated for $T_1 = 10\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ and no additional dephasing, while the dashed lines correspond to $T_1 = 17\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ and pure dephasing with $T_\phi = 16\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ due to low frequency noise (see text).](010206H1linesv3.eps){width="3.0in"}

As with all of the simulations discussed in this section so far, the solid lines in Fig. \[F010206H1lines\] were calculated with $T_1 = 17\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ and $T_\phi = 16\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$. This gives $T_2 = 10.9\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ and a Rabi decay time roughly equal to the two-level value[@Torrey49a; @Smith78a] of $T^\prime = {\ensuremath{\left[ 1/{\ensuremath{\left( 2 T_1 \right)}} + 1/{\ensuremath{\left( 2 T_2 \right)}} \right]}}^{-1} = 13.3\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$, consistent with measurements of the decay envelope of the escape rate. If pure dephasing were not present then $T^\prime$ would be 22.7 ns, which is significantly longer than what is observed. Although unlikely, our thermal measurement of $T_1 = 17\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ could be incorrect. To examine this possibility, simulations with $T_1 = 10\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ and $T_2 = 2 T_1$ (which also give $T^\prime = 13.3\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$) are shown with dotted lines in Fig. \[F010206H1lines\]. They are nearly identical to the curves calculated with dephasing, although [$\Gamma_{\infty}$]{} is somewhat smaller. Thus this data set alone cannot rule out dissipation-limited decoherence. However, the shorter $T_1$ reduces the prediction for [$\Gamma_{\infty}$]{} by roughly 15% over the full range of measured powers in Fig. \[F031806DGeq\], suggesting that additional dephasing is instead affecting the Rabi oscillations.

While dephasing is needed to faithfully reproduce features of the experimental measurements, its origin is unclear. In the simulations, we assume that each off-diagonal term of the density matrix decays with a common dephasing time of $T_\phi = 16\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ (along with decoherence due to dissipation, which is transition dependent). If this dephasing at high power were due to low frequency noise (similar to the effects of inhomogeneous broadening discussed before),[@Martinis03a; @Xu05b] the corresponding dephasing operator $D$ in [Eq. [(\[eqh\])]{}]{} could be expressed with $L = \sum {\ensuremath{\omega_{0 n}}} / {\ensuremath{\omega_{0 1}}} {\ensuremath{\left| n \right>}} {\ensuremath{\left< n \right|}}$ and $\lambda = 2 / T_\phi$ (similar to the harmonic oscillator number operator). With this choice of operator $L$ and strength $\lambda$, the [$0 \rightarrow 1$]{} dephasing is unchanged. Simulations with this sort of damping are shown with dashed lines in Fig. \[F010206H1lines\]. While the [$0 \rightarrow 1$]{} Rabi oscillations are modeled well, there is far too much decoherence for the higher order transitions. This should also be true for other implementations of low frequency noise that have dephasing rates that scale as $(d {\ensuremath{\omega_{n m}}}/d I_b)^2$. Thus having the ability to measure a wide range of transitions can reveal additional information about decoherence.

Discussion and Summary {#SSummary}
======================

We have presented measurements taken on a dc SQUID operated in such a way that one of the junctions behaves much like a simple current-biased junction. We find that the simplified Hamiltonian of [Eq. [(\[eqa\])]{}]{} gives an accurate description of the qubit dynamics. Our observation of subtle features of this model provides further confidence that it can be applied to design the gates needed for quantum computation. The analysis presented here can also be extended to describe the behavior of the other types of superconducting qubits mentioned in Sec. \[SIntro\].

We performed several checks on the model, including measuring the Rabi oscillation frequency for one and two-photon transitions and comparing the data with predictions based on the rotating wave Hamiltonian of [Eq. [(\[eqd\])]{}]{}. We find good agreement, despite the model not containing any damping. Thus the resonance shifts in multilevel Rabi oscillation frequencies are mainly determined by the anharmonic level structure and not the $T_1$ and $T_2$ of the qubit. These multilevel effects were identified under strong driving where it is clear that they could lead to errors in simple two-state rotations. Future experiments could mitigate these effects by proper pulse shaping and operating at low microwave power.[@Steffen03a; @Amin06a; @Lucero08a]

We also used a density-matrix simulation, which included the effects of decoherence, to calculate the tunneling escape rate. The results of the simulation agreed well with experimental Rabi oscillation data which were acquired with the device biased so that even the ground-state escape rate [$\Gamma_{0}$]{} was measurable. While this simple measurement (which ends with the qubit in the finite voltage state on every repetition) is not an ideal projective measurement to the qubit basis states [$\left| 0 \right>$]{} and [$\left| 1 \right>$]{}, it is particularly well-suited for the purposes of the current work. For one, it is extremely sensitive to excited state population due to a microwave drive or thermal transitions. The leakage information contained in Fig. \[F031806DGeq\] would be considerably more difficult to obtain from the pulsed single-shot measurements[@Cottet02a; @Claudon04a] that excel at determining the total population not in the ground state; nonetheless, these techniques have been successfully used to measure the second excited state population precisely.[@Palomaki07a; @Claudon07a; @Lucero08a] With the simple tunneling measurement, the large ratio ${\ensuremath{\Gamma_{n+1}}} / {\ensuremath{\Gamma_{n}}}$ provides a natural way of distinguishing which excited state is populated. In addition, as nothing is done to the qubit to initiate a measurement, this method is not limited by a measurement fidelity and does not suffer from the potential problems associated with changing the bias location in order to perform state readout.[@Cooper04a; @Zhang06a; @Palomaki07a; @Claudon07a]

The main drawback to the escape rate measurement is that it does not directly produce individual state populations. However, the density-matrix simulations do provide each level’s contribution to the total escape rate and its occupation probability. An example for a [$0 \rightarrow 1$]{} oscillation is shown in Figs. \[F032106FN10\](b) and \[F032106FN10\](c). The plots suggest that for high power, the escape rate during a Rabi oscillation is dominated by the contributions $\rho_{22} {\ensuremath{\Gamma_{2}}}$ and $\rho_{33} {\ensuremath{\Gamma_{3}}}$ from states [$\left| 2 \right>$]{} and [$\left| 3 \right>$]{}. The figure also shows that $\rho_{22} {\ensuremath{\left( t \right)}}$ and $\rho_{33} {\ensuremath{\left( t \right)}}$ have nearly the same form as $\rho_{11} {\ensuremath{\left( t \right)}}$, albeit with a much smaller oscillation amplitude.[@Amin06a] Thus changes in $\Gamma$ are reflective of the underlying oscillation of [$\left| 1 \right>$]{}, which is why we believe the frequency analysis in Fig. \[F031806Dfosc\] is valid.

The populations in Fig. \[F032106FN10\](c) have been normalized at each time $t$. While the probability of the junction leaving the supercurrent state during a microwave pulse depends strongly on [$I_{r\!f}$]{} and the bias conditions, $\rho_{tot} = 0.42$ at $t = 30\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ for this data set. Thus, while we wish to use tunneling as a probe of $\rho_{nn}$, it is clearly playing a large role in the evolution of the system. However, repeating the density-matrix simulation in the absence of tunneling (*i.e.*, $G {\ensuremath{\rho}}= 0$) results in a less than 5% change in the normalized populations of [$\left| 0 \right>$]{} and [$\left| 1 \right>$]{} during the oscillation. The subsequent decay is significantly different without the fast decays due to escape. These insights from the simulation hinge on the assumption that tunneling does not affect $T_1$ and $T_\phi$.

A few remarks should be made about the many input parameters required for the simulations. For example, the seven-level density-matrix calculation shown in Fig. \[F010206H1\](c) needed seven escape rates, six energy-level spacings, and 28 matrix elements of $\gamma_1$, all of which are functions of $I_b$. We determined these from the properties of a tilted washboard potential defined by $I_{01}$ and $C_1$, which were measured independently by low power spectroscopy. To accurately reproduce the measured escape rate (particularly at low power), we also had to add several other features to the simulation, each motivated by separate measurements. In particular, we added a microwave source at [$\omega_{0 2}$]{} to mimic noise (whose magnitude was found from the microwave-free escape rate), inhomogeneous broadening due to bias noise (estimated from spectroscopic resonance widths), and a finite time resolution (consistent with the bandwidth of the detection electronics). The dissipation time $T_1 = 17\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ was estimated by an independent thermal escape rate experiment,[@Dutta04a] while the dephasing time $T_\phi = 16\ {\ensuremath{{\ensuremath{\mathrm{ns}}}}}$ (and the way it was incorporated into the simulation) was chosen to maximize the agreement with the Rabi data. For simplicity, both time constants were assumed to be independent of frequency. Removing this condition in the analysis could yield additional information about the decoherence in the system, as the oscillations are most sensitive to noise at the Rabi frequency.[@Martinis03a; @Ithier05a] Finally, the conversion between the microwave power at its generator and the current through the qubit junction was itself calibrated by the observed oscillation frequency, and the resulting factor was allowed to be a function of microwave frequency. Thus a wide range of multilevel phenomena was explained with just a few truly free parameters.

In addition, we found that the output of the simulations (whether that be an oscillation frequency or escape rate) converged as the number of levels $N$ was increased. For any of the five or seven-level simulations discussed here, an additional level produced a change too small to be detected by the experiment; more levels were required at high power or low $I_b$ for the solution to converge satisfactorily. Surprisingly, the highest included levels are predicted to lie above the barrier. However, small errors in the energy levels or matrix elements for the highest levels do not produce large changes in the final escape rates.

Certain features of the data were not reproduced by the simulations. Although not discussed here, the transition spectra of our Nb qubits show a series of splittings (all less than 10 MHz wide), similar to but smaller than those reported in other superconducting qubits.[@Simmonds04a; @Plourde05a; @Claudon07a; @Dong07a] The extra degrees of freedom responsible for these features were not included in the device Hamiltonian, except perhaps in some effective way through $T_1$ and $T_\phi$ of the density-matrix simulation. For weak coupling and high microwave power, the impact of individual two-level systems on Rabi oscillations is expected to be small.[@Meier05a; @Ashhab06a] However, it is unclear if a bath of quantum systems could be responsible for the anomalous decay in Fig. \[F032106FN10\]. As other groups have observed strong dissipation in Nb trilayer junctions,[@Martinis02a; @Lisenfeld07a] it is likely that the wiring insulation or other details of the fabrication are playing a critical role.[@Martinis05a] While further work is needed to identify (and potentially eliminate) the actual microscopic sources of decoherence in this system, the multilevel features in the dynamics of our qubit are largely understood.

We have benefited greatly from discussions with B. K. Cooper, P. R. Johnson, H. Kwon, J. Matthews, K. D. Osborn, B. S. Palmer, A. J. Przybysz, R. C. Ramos, C. P. Vlahacos, and H. Xu. The work was funded by the National Science Foundation through the QuBIC Program, the National Security Agency, and the state of Maryland through the Center for Nanophysics and Advanced Materials, formerly the Center for Superconductivity Research.

[66]{} natexlab\#1[\#1]{}bibnamefont \#1[\#1]{}bibfnamefont \#1[\#1]{}citenamefont \#1[\#1]{}url \#1[`#1`]{}urlprefix\[2\][\#2]{} \[2\]\[\][[\#2](#2)]{}

, , , ****, ().

, , , .

, ****(11), ().

, ****, () .

, ****, ().

, , , , , , ****, ().

, , , , , ****, ().

, , , , , , ****, ().

, , , , ****, ().

, , , , ****, ().

, , , , , , , ****, ().

, , , , , , , , , , , ****, ().

, , , , , ****, ().

, , , , , , , , , , , ****, ().

, , , , , , , , , ****, ().

, , , , , , , , , , ****, ().

, , , , ****, ().

, , , ****, ().

, , , , ****, ().

, , , , , , ****, ().

, , , ****, ().

, , , , , , ****, ().

, , , , , ****, ().

, , , ****, ().

, , , , , ****, ().

, ****, ().

, , , ****, ().

, ****, ().

, ****, () .

, , , , , ****, ().

, , , , ****, ().

, , , , , , , , , , ****, ().

, , , , , , , , , , , ****, ().

, , , , , , ****, ().

, ****, ().

, in **, edited by , , (, , ).

, ****, () .

, , , ****, ().

, ****, ().

, Ph.D. thesis, , .

, ****, ().

, ****, ().

, , , , , , , , , , , ****, ().

, , , ****, ().

, , , , , , , ****, ().

, , , , , , , , ****, ().

, , , , .

, , , , , , , , ****, (); .

, , , , ****, ().

, ****, ().

, ** (, , ).

, , , , ****, ().

, , , ****, ().

, , , , , , , , , , ****, ().

, , , , , , , ****, ().

, ****, ().

, , , , , ****, ().

, , , , , , ****, ().

, , , , , , , , ****, ().

, , , , ****, ().

, , , , , , , , , , , ****, ().

, , , , , , , ****, ().

, , , , , , , , , , ****, ().

, , , ****, ().

, , , , , , , , , , , ****, ().
