Theoretical study of reflection spectroscopy for superconducting quantum parametrons

Superconducting parametrons in the single-photon Kerr regime, also called KPOs, have been attracting increasing attention in terms of their applications to quantum annealing and universal quantum computation. It is of practical importance to obtain information of superconducting parametrons operating under an oscillating pump field. Spectroscopy can provide information of a superconducting parametron under examination, such as energy level structure, and also useful information for calibration of the pump field. We theoretically study the reflection spectroscopy of superconducting parametrons, and develop a method to obtain the reflection coefficient. We present formulae of the reflection coefficient, the nominal external and the internal decay rates, and examine the obtained spectra. It is shown that the difference of the populations of energy levels manifests itself as a dip or peak in the amplitude of the reflection coefficient, and one can directly extract the coupling strength between the energy levels by measuring the nominal decay rates when the pump field is sufficiently large.


Introduction
Classical parametric phase-locked oscillators [1], called parametrons [2], were operated as classical bits in digital computers in the 1950s and 1960s. Recently, parametrons in the single-photon Kerr regime [3], where the nonlinearity is larger than the decay rate, have been attracting much attention in terms of their applications to quantum information processing. Parametrons were applied to the qubit readout [4,5] in the circuit QED architecture, a promising platform of quantum information processing [6,7,8,9,10,11,12,13]. Quantum annealing [14,15,16,17] and universal quantum computation [18], which utilize the quantum nature of parametrons in a superconducting circuit, have been proposed. More recently, bias-preserving gates [19] and singlequbit operations [20] were studied theoretically and experimentally, and the exponential increase of the bit-flip time with the cat size was observed [21].
A parametron in the single-photon Kerr regime is operated by an oscillating pump field. Therefore, to obtain information of the parametron under the pump field and to calibrate the amplitude of the pump field are practically important. In previous studies, state tomography of parametrons using the power spectrum density were demonstrated [22]; energy differences between either of the two highest energy levels and a lower energy level of a parametron were measured by mapping the parametron to a Fock qubit [20]; and microwave responses of parametrons without a pump field were experimentally investigated in a wide range of the Kerr nonlinearity [23]. However, theories of reflection and transmission spectroscopy of parametrons have not been developed in spite of the fact that they are important and routinely applied to resonators to examine the energy level structure and their quality. Spectroscopy will provide useful information of superconducting parametrons under examination, such as energy level structure.
Recently, a fast and accurate gate operation of a parametron utilizing energy levels outside of the qubit space was proposed [24]. In the method, the couplings between either of the two highest levels and other lower levels induced by a drive pulse are essential. For implementation of the technique, it is important to measure the energy differences and the couplings between relevant levels.
In this paper, we develop a method to obtain the reflection coefficient of a pumped superconducting parametron. We present simple formulae of the reflection coefficient and the nominal decay rates ‡. It is shown that one can directly extract the amplitude of the coupling coefficients between energy levels by measuring the nominal decay rates when the pump field is sufficiently large. Moreover, we show that the nominal internal decay rate increases with the pump strength and eventually exceeds the nominal external decay rate even if the original internal decay rate of the parametron without a pump field is negligible. Our method of spectroscopy does not require pulsed operations of a parametron, and can be implemented with a standard reflection-measurement setup routinely used for superconducting circuit QED architectures.

Model
We develop a theory to obtain the reflection coefficient for parametrons. Our method is similar to that in Ref. [25] developed for a driven three-level system. We consider a system composed of a parametron attached to a transmission line (TL) as depicted in figure 1(a). An incoming and an outgoing microwaves propagate in the same TL. The ‡ The external and the internal decay rates or the quality factors of resonators are routinely measured by spectroscopic methods. They are obtained by fitting the spectra to an analytic form. The nominal decay rates of a parametron can be obtained in the same manner as a resonator, and can provide information of the measured parametron as explained in Sec 4. parametron is pumped by an external oscillating magnetic flux Φ(t). An external oscillating magnetic flux Φ(t) is used to pump the parametron. ω is the resonance frequency of the parametron when no pump filed is applied. The parametron is located at r = 0. (b) Effective model of (a). In the theory in section 3, the negative and positive regions are assigned to the incoming and outgoing fields, respectively.

Hamiltonian of the system is given by
where the first line is Hamiltonian of the parametron under a pump field [22] (see Appendix A for derivation), and the second line is Hamiltonian of the eigenmodes of the TL and the coupling between the modes and the parametron. The third line is Hamiltonian of the eigenmodes representing a loss channel and their coupling to the parametron. The decay to the channel represents the internal decay of the parametron. The second term in the right-hand side of equation (1) gives rise to an anharmonic term in a rotating frame as shown later. The third term originates from the pump field Φ(t). Here, β and ω p are the amplitude and the angular frequency of the pump field, respectively. The annihilation operator of the parametron and those of the eigenmodes of the TL (loss channel) with the wave number k are denoted by a and b k (c k ), respectively. The decay rate to the TL (loss channel) and the phase velocity in the TL (loss channel) are denoted by κ ex (κ int ) and v b (v c ), respectively. Hereafter, we assume v c = v b .

Reflection coefficient
In this section, we show the method to calculate the reflection coefficient of a parametron. In our method, we use the input-output relation. We derive the inputoutput relation to make this paper self-contained, although the derivation is based on a standard approach and can be found, e.g., in Ref. [5]. For a parametron, energy eigenstates in a rotating frame are important because two of them are utilized as qubit states. We develop a method to obtain the reflection coefficient for a parametron using its energy eigenstates as a basis in section 3.1. This method enables one to obtain the reflection coefficient without integrating master equation. The Heisenberg equation of motion of b k is represented as A formal solution of equation (2) is where t > 0. We formally extend the lower limit of k to −∞ in order to introduce the real-space represention of the field operator defined bỹ where r runs over −∞ < r < ∞. The negative and positive regions are assigned to the incoming and outgoing fields, respectively, as depicted in figure 1(b). Thus, the input field operatorb (in) r and the output field operatorb The introduction of the real-space representation has been validated in reference [26]. Using equations (3) and (4), we can obtaiñ where θ is the Heaviside step function. Using equation (5) with r = 0 in equation (6), we obtain the input-output relation, where we used θ(0) = 1/2. The Heisenberg equation of motion of a with equations (4), (5) and (7) leads to where a abbreviates a(t), and we formally extended the lower limit of k to −∞ in equation (1). Here,c  v b t , and κ tot = κ ex + κ int . Now, we assume that an input microwave is applied from the TL attached to the parametron. We consider a continuous mode version of a coherent state: with the overall vacuum state |v and a normalization factor N . Considering that the input wave propagates in the positive-r direction, E in (r) represents the input microwave at the initial moment as given by where E and ω in are the amplitude and the angular frequency of the incoming microwave, respectively. We assume that at the initial time the parametron is unexcited, and the input microwave has not arrived at the parametron yet. Because |Ψ i is a coherent state, it is an eigenstate of the initial field operatorb r (0). We can obtaiñ using equations (9) and (10). We multiply equation (6) by e iωpt/2 with r = +0 and take the expectation value with respect to |Ψ i to obtain whereB (out) and Equation (11) leads to In this paper, we focus on the Fourier component of B (out) +0 (t) with a frequency of −ω in + ω p /2, which is the same as the frequency of the input field, although our formalism can be used to obtain other frequency components of the reflected field. We define the reflection coefficient as where B (out) (16) can be rewritten using equations (12) and (15) as where

Representation of reflection coefficient with density matrix elements
We consider equations of motion of the system under consideration to calculate the reflection coefficient. The time evolution of A is governed by where we used equations (8) and (11) and omitted rapidly oscillating terms (rotating wave approximation). The approximation is valid when ω p ≫ β, ∆, χ. Effects of these rapidly oscillating terms on controls of a parametron were studied in Ref. [27]. Here, ∆ is the detuning defined by The master equation, which leads to the same equations of motion, is represented as with where ρ is the density operator.
In order to derive an analytic formula of the reflection coefficient, we use energy eigenstates |φ m of H 0 in equation (22) to rewrite A as where ρ nm = φ n |ρ|φ m , and X mn = φ m |A|φ n . Using equations (17) and (24), the reflection coefficient can be represented as where ρ (25) is the contribution to the reflection coefficient from the transition from |φ m to |φ n .
The equation of motion of the density matrix element is written aṡ where ω m is an eigenvalue of H 0 / ; Y mn = φ m |A † A|φ n ; and Ω is defined by The Fourier transform of equation (26) with a frequency of −ω in + ω p /2 (= −ω in ) leads to These equations are used to obtain the density matrix elements in equation (25) as shown in the following section.

Weak input field limit
In principle, the reflection coefficient in equation (17) can be obtained using the density matrix which can be calculated by integrating the master equation (20). However, it is time consuming to integrate the master equation for sufficiently long time for a parametron. In this section, we present an alternative method providing an approximate reflection coefficient in the weak input field regime, where the diagonal elements of the density matrix ρ are approximately the same as those for the stationary state without input field. The method does not require to integrate the master equation. The effect of the finite input field can be also taken into account, which will be discussed in Appendix B.
We consider the case that the input field is nearly resonant with the transition from |φ m to |φ n , that is, ω in − ω p /2 + ω m − ω n ≃ 0. We assume that the off-diagonal element of the density matrix ρ We also assume that the diagonal elements are the same as the stationary state without input field because the input field is sufficiently weak. Then, ρ (F) nm [−ω in ] can be obtained using equation (28) as where ∆ nm = ω in − ω p /2 − ω n + ω m . We used X mm = 0 in equation (29), which is valid because |φ m has even or odd parity, that is, |φ m is a superposition of even-photonnumber states or odd-photon-number states. Using equation (29) in equation (25), the reflection coefficient is approximately represented by the following simple form, Note that the diagonal elements of the density matrix in equation (31) are for the stationary state of the case without input field. When either of the resonant energy levels is occupied, the amplitude of the reflection coefficient changes from unity, and thus a dip or peak in the reflection coefficient appears.
Around a dip or peak corresponding to the transition from |φ m to |φ n , Γ is approximated by Γ mn defined by On the other hand, the reflection coefficient of a linear resonator is represented as [29] Γ r = 1 + κ where κ can be interpreted as the nominal external and internal decay rates, respectively.
The internal and the external decay rates of a linear resonator can be obtained via a fitting of the measured reflection coefficient to the analytic form of equation (34). The nominal decay rates of the parametrons can be obtained in the same manner with equation (33), and can provide information of the measured parametron such as ) and Y mm + Y nn . Note that κ ex and κ int in equation (33) can be obtained via measurements of the parametron without pump field.
The measurement ofκ mn ex is rather useful in the strong pump regime, where β is sufficiently larger than χ and ∆. It is known that the stationary state is the maximally mixed state of the two highest levels in the strong pump regime [14]. For example, in the case of ∆ ≤ 0, we have ρ  (33) gives the amplitude of the coupling coefficient |X mn | between either of the two highest levels and a lower level as |X mn | = 2κ mn ex /κ ex , where m = 0, 1 and n ≥ 2. Therefore, the amplitude of the coupling coefficients can be directly extracted by the measurement ofκ mn ex . Recently, a fast gate operation of a parametron utilizing energy levels outside of the qubit space was proposed [24]. It is useful to experimentally extract |X mn | for tailoring a control field in such protocols to improve the gate fidelity.
The reflection coefficient can be calculated also by a straightforward but time consuming manner of integrating the master equation (20). In Appendix C, we compare the results of the two methods to numerically verify the approximate method.  (31) were numerically calculated using equation (26). The used parameter set is: χ/2π = 30 MHz, κ ex /2π = 0.4 MHz and κ int /2π = 4 MHz, and is typical for superconducting parametrons [22,23]. When β ≃ 0, there is only a single dip of |Γ| at ω in − ω p /2 = ∆ corresponding to the transition from Fock state |0 to Fock state |1 as seen in figure 2. On the other hand, the spectra show the interesting behaviors as β increases: the dip corresponding to the transition |0 → |1 disappears while other peaks and dips appear; the frequencies corresponding to some of the peaks and dips increase with β, while the others decrease; a dip (peak) changes to a peak (dip) in figure 2(d). And, the pattern of the spectrum depends on the detuning. In the following, these behaviors of |Γ| are explained. Each dip and peak of |Γ| corresponds to a transition between eigenstates of H 0 in equation (22). Therefore, it is useful to examine the eigenvalues of H 0 shown in figure 3. Each eigenstate of H 0 denoted by |m coincides with Fock state |m when β = 0. We denote the eigenenergy as ωm. Figure 4 shows the same spectra as figure 2 with the curves representing energy differences, ∆ωñm = ωñ−ωm, corresponding to the transition, |m → |ñ . The dips and peaks match to the curves for the energy differences. Thus, the spectra reflect the energy level structure of H 0 . Therefore, the spectra can provide information on the energy level structure of the parametron. As seen in figure 3, the order of levels, |m , depends on the detuning [28]. Thus, the pattern of the spectrum also changes depending on the detuning. For example, the order of the distinct dip and peak around ω in = ω p /2 are opposite in figures 4(a) and 4(c). It reflects the difference in the order of |0 and |1 observed in figures 3(a) and 3(c). As seen from equation (31), the finite population difference (|ρ
It is useful to examine ηmñ = −|Xmñ| 2 (ρ (F) (31) to explain the appearance and disappearance of dips and peaks of |Γ|. Note that Ymm and Yññ are always positive, and ∆ñm is zero at the resonance in equation (31). If ηmñ is positive, a peak corresponding to the transition |m → |ñ shows up. If negative, a dip appears. Figure 7 shows ηmñ as a function of β for four sets of (m,ñ). The results for ∆/2π = −7, 0, 7 MHz indicate: 1. the dip for the transition |0 → |1 vanishes as β increases; 2. the dip for |1 → |2 and the peak for |1 → |0 appear as β increases, although the peak vanishes when β is further increased; 3. transition |2 → |1 is hardly seen. This is because that the coupling coefficient, X21, is approximately zero as shown in figures 6(a)−(c). The results for ∆/2π = 20 MHz indicate: 1. the dip for |1 → |2 and the peak for |2 → |1 appear for intermediate value of β; 2. the transition |0 → |1 gives rise to a dip of |Γ| for relatively small β and a peak for β/2π ≃ 10 MHz, and then vanishes as β is further increased; 3. the transition |1 → |0 gives rise to a peak for relatively small β and a dip for relatively large β. The change of the sign of ηmñ corresponding to the transitons |1 ↔ |0 comes from the crossing of the populations of |0 and |1 observed in figure 5(d). Thus, these results explain the profile of the  spectrum in figure 2. Figure 8 shows the nominal external and the nominal internal decay rates in equation (33). The nominal external decay rate,κ (01) ex , decreases to zero as β increases. Some of other nominal external decay rates become finite for β = 0 although they are approximately zero for β ≃ 0. On the other hand, nominal internal decay rate increases rapidly with respect to β. Thus, the dips and peaks tend to broaden as β increases. Even if the original internal decay rate of the parametron without a pump field is negligible, the nominal internal decay rate increases with the pump strength and eventually exceeds the nominal external decay rate (see Appendix D).
The Wigner function in figure 9 illustrates the stationary state for the values of β indicated by the triangles in figure 3(a). The profile of the Wigner function depends on the detuning for relatively small value of β although it is insensitive to the detuning for β ≫ |∆|, χ as seen in figures 9(d), (h), (l), (p). This is because the stationary state becomes the maximally mixed state of the coherent states, |α and | − α , with α ≃ 2β/χ for β ≫ |∆|, χ.

Summary
We have theoretically studied the reflection spectroscopy of a pumped superconducting parametron. We have developed a method to obtain the reflection coefficient of a parametron and have derived formulae of the reflection coefficient, the nominal external and internal decay rates. This method can also take into account the effect of the input field beyond the limit of weak input field. It has been shown that the peak or dip can appear in the amplitude of the reflection coefficient when there is finite difference between the populations of energy levels resonantly coupled by an input field, and the sign of the difference determines whether we have a dip or peak. We have shown that the nominal internal decay rate increases with the pump strength and eventually exceeds the nominal external decay rate even if the original internal decay rate of the parametron without a pump field is negligible. The obtained spectrum provides information of the superconducting parametron, such as energy level structure and amplitude of coupling coefficients between energy levels, and also useful information for calibration of the pump field.

Appendix A. Hamiltonian of a parametron
We derive an effective Hamiltonian for a parametron to make this paper self-contained although it can be found in Ref. [22]. We consider a parametron consisting of a SQUIDarray resonator with N SQUIDs depicted in figure A1. The effective Hamiltonian of the system is given by where φ and n are the overall phase across the junction array and its conjugate variable, respectively. E J is the Josephson energy of a SQUID. We assume that all the Josephson junctions are identical. The effective Hamiltonian (A.1) with a single degree of freedom, φ, is valid provided that E J is much larger than the charging energy of a single junction [31,32]. E C is the charging energy of the resonator, including the contributions of the junction capacitances C J and the shunt capacitance C, and can be experimentally extracted or calculated by finite-element capacitance simulation [22]. The Josephson energy is modulated as E J (t) = E J + δE J cos ω p t by the external magnetic flux, Φ(t), threading the SQUIDs. We take into account up to the fourth order of φ/N in equation (A.1) to obtain an approximate Hamiltonian where ω = 1 8E C E J /N, χ = E C / N 2 and β = ωδE J /8E J . Here, β is called amplitude of the pump field in the main text. n and φ are related to the annihilation operator a as Above, we considered the parameter regime, where φ 0 /N = 2 χ/ω is sufficiently smaller than unity so that the expansion of cos(φ/N) is valid, and took into account up to the fourth order of φ/N to see the effect of the Kerr nonlinearity. In equation (A.2), we neglect the last term assuming that χβ ≪ ω, and drop c-valued terms to obtain the following Hamiltonian H = ωa † a − χ 12 (a + a † ) 4 + 2β(a + a † ) 2 cos ω p t.

Appendix B. Effect of input field
In the main text, we considered the weak input field limit. The diagonal elements of the density matrix, ρ (F) mm [0], were calculated assuming that they are not changed by the input field, Ω = √ v b κ ex E. However, we can take into account the effect of finite Ω by solving the Fourier transform of equation (26) to obtain the elements of the density matrix. In the numerical simulations of this section, we assume ρ  Figure B1 shows the amplitude of the reflection coefficient in equation (25) as a function of ω in and β. The result for Ω/2π =1 MHz is approximately the same as the results in figure 2 for the weak input field limit. The peaks (dips) become low (shallow) for larger Ω. This is attributed to thatκ

Appendix C. Direct numerical simulations with integration of master equation
In section 4, an approximate formula for reflection coefficient was derived. The reflection coefficient can be calculated also by a straightforward but time consuming manner. We integrate the master equation (20) to obtain the density matrix and calculate A [−ω in + ω p /2]. Equation (17) is used to obtain the reflection coefficient. We compare the results with those obtained by the method in section 4.
In the master equation, we set the initial state of the parametron to the stationary state. We integrate the master equation for 0 ≤ t ≤ 440 ns with the constant input field. We used a fourth-order Runge-Kutta integrator with the time step of less than 0.012 ps. Figure C1 shows the amplitude of the reflection coefficient for the both methods. The dip at (ω in − ω p )/2π ≃ −53 MHz (-94 MHz) corresponds to the transition |1 → |2 (|0 → |3 ). It is seen that the results for the method in section 4 approximates well especially near the resonance (dips). There is a discrepancy between the two results between the two dips (see figure C1(a)), which we attribute to the fact that we neglect the interference between different transitions in the method in section 4. This discrepancy becomes small when we decreases κ int as the dips are well separated as seen in figure C1(b,c). The nominal internal decay rate increases with respect to the pump amplitude even if the original internal decay rate of the parametron without a pump field is negligible as shown in this section. To observe this fact, we consider a fictitious case where κ int is zero. Figure D1 shows the amplitude of the reflection coefficient calculated for the weak input field limit as a function of ω in and β for κ int = 0. The dips and peaks are sharper than those in figure 2 for κ int /2π = 4 MHz. Figure D2 shows the nominal external and the nominal internal decay rates in equation (33). The nominal external decay rate is approximately the same as that in figure 8 for κ int /2π = 4 MHz. We attribute this to the fact that the diagonal elements of the density matrix are approximately the same in both cases. On the other hand, nominal internal decay rate is much smaller than that in figure 8 due to vanishing κ int . However, the nominal internal decay rate increases rapidly with respect to β because Ymm and Yññ increase with β. Thus, the broadening of the dips and peaks occur also in the case of zero κ int .