Mediated tunable coupling of flux qubits

It is sketched how a monostable rf- or dc-SQUID can mediate an inductive coupling between two adjacent flux qubits. The nontrivial dependence of the SQUID's susceptibility on external flux makes it possible to continuously tune the induced coupling from antiferromagnetic (AF) to ferromagnetic (FM). In particular, for suitable parameters, the induced FM coupling can be sufficiently large to overcome any possible direct AF inductive coupling between the qubits. The main features follow from a classical analysis of the multi-qubit potential. A fully quantum treatment yields similar results, but with a modified expression for the SQUID susceptibility. Since the latter is exact, it can also be used to evaluate the susceptibility--or, equivalently, energy-level curvature--of an isolated rf-SQUID for larger shielding and at degenerate flux bias, i.e., a (bistable) qubit. The result is compared to the standard two-level (pseudospin) treatment of the anticrossing, and the ensuing conclusions are verified numerically.


I. INTRODUCTION
Superconducting devices made of mesoscopic Josephson junctions are one of the few remaining viable candidates for the realization of quantum bits (qubits). The major types encode quantum information in the charge, phase, or flux degree of freedom [1]. Quantum coherence in such devices was conclusively demonstrated at the turn of the century. While much remains to be done in improving the coherence times, control, and readout of individual qubits, another key challenge on the road towards quantum-register prototypes is the coupling of several qubits. In the past years, coupling was achieved for each of the mentioned qubit types [2,3].
In all published experiments, the coupling strength is a constant determined at the time of fabrication. For each model of quantum computing (gate-operations, adiabatic/groundstate, etc.), however, it would be very desirable to have the coupling tunable in situ, even if this requirement can sometimes be circumvented in "software" [4]. In particular, whenever the coupling can be continuously tuned between two values with opposite signs, it can be switched off naturally.
For charge qubits, a design fulfilling these requirements was presented in Ref. [5]. The coupling between the two qubits is mediated via a third device, the bias parameters of which can be used to tune the coupling constant. While the coupling element thus is nontrivial, it is important that it be passive, i.e., that it does not get entangled with the qubits so that the system's total effective Hilbert space remains four-dimensional. In practice, this means that the coupler must stay in its ground state, adiabatically following the qubit dynamics.
Mediated coupling can be adapted to flux qubits [6]. This paper's purpose is to contribute to its analytical theory. The main results can be obtained from a purely classical analysis of the system's potential, which may have some independent pedagogical interest. This is carried out in Sec. II for two choices of the coupler, both SQUID devices. A full quantum analysis is presented in Sec. III; the main result for the tunable coupling constant is highly analogous to its classical counterpart, but some additional effects are uncovered concerning a renormalization of the qubit parameters. In the formula for the coupling constant, the nontrivial factor is an exact expression for the coupler's ground-state magnetic susceptibility. In Sec. IV, this is applied to a different regime, where the coupler behaves like a qubit itself, allowing for a comparison with the predictions of the two-state model. Some concluding remarks are made in Sec. V.

A. rf-SQUID coupler
Consider three rf-SQUIDs a through c in a row (Fig. 1). The a-and c-devices are supposed to work as (effectively) degenerately biased flux qubits, while the b-SQUID is a coupling element, tunable by an external flux bias. To derive the indirect a-c coupling it suffices to consider the potential, consisting of both Josephson and magnetic terms ( = 1 throughout), FIG. 1: A system of three rf-SQUIDs, of which the two outer ones function as qubits and the inner one as a coupler.
with the inductance matrix That is, the direct AF inductive a-c coupling is ignored. In principle, this can always be achieved through a gradiometric layout. However, as long as M ac is small, one can probably just as well simply add an interaction M ac I a I c to the final result (9) below. In accordance with the roles played by the various loops, the bias fluxes (in phase units) are chosen as Thus, φ x a,c compensate externally for the shielding flux which the b-loop couples into the a-and c-SQUIDs; see (6) for the precise definition of φ The b-SQUID should act as a passive coupler without introducing additional degrees of freedom, so bistability must be avoided. This can be achieved either by biasing it with a flux close to an integer number of quanta (φ x b small), or for any flux bias by taking the shielding . The calculation that follows is valid in either case, and their relative merits will be discussed afterwards. 1 We proceed by expansion in M. Without suggesting that the regime of M/L ≪ 1 (say, distant loops) is the most practical, this leads to a transparent result, which can guide numerical studies in the general case. The b-mediated a-c coupling is expected to be O(M 2 ), so the junction phases will be written as these can be determined by solving the equilibrium condition ∇ φ U = 0. In leading order, the φ (0) j are just the stationary phases of an isolated rf-SQUID, where φ (0) a,c − π can have either sign. To the same order as (4), one has (7) While (7) may look tedious, using it consistently actually leads to significant cancellation below.
In O(M), one finds that φ (1) a,c = 0 due to the special choice of φ x in (3), while the b-loop picks up the shielding fluxes of the neighboring ones, with a prefactor reflecting its susceptibility: Calculating the φ is straightforward but unnecessary: since one expands around a minimum of U, the φ (2) j do not contribute to the relevant order. All that remains is to substitute (3)-(8) into U in (1). Since, e.g., (φ (0) a −π) 2 does not depend on the qubit state, one is left with The product of mutual inductances is expected geometrically, while I a I c ∝ σ z a σ z c in qubit parlance. The fraction is worth a closer look: Thus, for small φ x b , the coupling is AF, but it changes sign to FM as φ x b → π (only attainable for β b < 1).
This behaviour, and the coupling mechanism in general, have a transparent interpretation. Consider the self-flux curve φ Fig. 2), which up to constants simply is the shielding current versus applied flux. If the a-qubit flips from ↑ to ↓, then due to the mutual inductance −M ab < 0 this effectively increases φ x b . For small φ x b , the curve has slope < 0, so that this increase generates a self-flux in the ↓-direction. In its turn, the latter, through −M bc < 0, acts as an ↑-bias for the c-loop, favouring the ↑-state there, i.e., opposite to the state of the a-qubit. This also explains the bound of unity exemplified by for an rf-SQUID coupler with β b = 0.95, according to Eq. (6).
the second line of (10), since the maximum AF response is perfect shielding for β b → ∞ (an uninterrupted superconducting loop). However, near φ x b = π the argument works in the opposite direction. There, the self-flux curve has slope > 0 so that, differentially, the self-flux does not shield at all but actually cooperates with the external flux increase. Thus, the coupling changes sign to FM. Moreover, as β b ↑ 1, the slope of this curve increases without bound (a precursor to bistability), so that the FM coupling can actually be large; this corresponds to the zero in the denominator of (10). On the one hand, the divergence is never realized in practice, since for any finite M's one deals with finite differences, not slopes, on this curve (a nice counterpart to the discussion in [5] for charge devices 2 ). On the other, this potentially large FM coupling makes it possible to overcome any residual direct AF coupling through M ac . Since (10) shows that the large-β b regime is quite inflexible, and that for β b ≪ 1 the coupling strength is limited by the small shielding flux, this case β b < ∼ 1 is thought to be the preferred one. The above is mostly a straightforward derivation from (1), so very few specifics of the circuit are involved. For instance, the a-and c-SQUIDs could also be placed inside a large b-loop, which changes the signs of both M ab and M bc , so that (9) is invariant. This design modification makes it clearer that the b-loop is mostly a flux transformer [7], with a weak link providing the tunability. Moreover, the final result (9) depends only trivially on the properties of the a-and c-devices-through the flux which these couple into the b-SQUID. Hence, the generalization to other types of flux qubits should be obvious.

B. dc-SQUID coupler
The above indicates that the simple rf-SQUID should be a satisfactory coupling element. Still, it is worthwhile to also examine a dc-SQUID [6], which provides additional flexibility since it can be both flux-and current-biased. In particular, if only a current bias I x b would turn out to suffice for coupling tunability, that could be preferable in situations where a flux bias is problematic due to stray fields etc. The b-device now has separate left and right arms, designated as b1 and b2, together carrying I x b (Fig. 3). Again ignoring the direct a-c mutual inductance (not that this makes much difference anywhere), we thus take the coil-to-coil inductance matrix as  However, always I b2 = I b1 + I x b , so the flux-current equation relates three-vectors pertaining to the loops: with On the other hand, the four junction phases are all independent dynamical variables. They are related to the fluxes by (φ a , φ b1 +φ b2 , φ c ) ⊤ = 2eΦ. One should first of all determine the Hamiltonian. Of the Kirchhoff (circuit) equations, the first four are simply the Josephson relationsφ a = 2eQ a /C a etc., with Q a the junction charge and C a its capacitance. The other four express current conservation: It is readily seen that these are the Hamilton equationsφ j = 2e∂H/∂Q j andQ j = −2e∂H/∂φ j , for with j ∈ {a, b1, b2, c}. It is less apparent that H is actually symmetric in the two SQUID arms, as it should be: H is a function of the junction phases and charges only, and the asymmetric choice of I b1 as the loop current was merely an interim convention. That is, designating the last two terms of (23) as "magnetic" and "bias" energy, respectively, is a bit arbitrary. Working out L −1 , one sees that up to a constant, one can rewrite It would be instructive to see if the above, especially the form (25) for H bias , can be reproduced (presumably almost mechanically) by the "node formalism" of circuit analysis [8,9].
From here on, the analysis is a straightforward generalization of the rf-SQUID case, again writing the phases as in (4), and expanding H bias in (25) to the same order. We again want the a-and c-qubits to be effectively at degeneracy, which presently means that their own external flux should compensate also for I x b : In O(M 0 ), the stationary phases obey the standard equations for the isolated devices: Here, β a(c) ≡ 4e 2 L a(c) E a(c) is standard, but note that the definitions β bj ≡ 4e 2 L b E bj (j = 1, 2) involve the full loop inductance of the b-SQUID, not the individual arm inductances.
In O(M), our choice (26) for φ x again ensures that φ c = 0, while the equations for φ (1) b1 (2) are coupled, as were (28)-(29) for φ (0) b1(2) (although, unlike the latter, they are of course linear): Like for the rf-SQUID, it is best to leave the φ (2) j unevaluated in the M-expansion of U, and one sees that their contribution cancels (to this order), since one expands around a potential minimum. Taking advantage of some further cancellation and dropping all terms which do not depend on the qubit state, one arrives at in which all tunability thus comes via the dependence of φ ) as given in (28)-(29). For a simple consistency check, let us take I x b → 0 and E b2 → ∞. Then, (29) shows that φ b1 ; in view of the remark below (29), this is exactly the result for rf-SQUID coupling.
Another well-known case is the symmetric dc-SQUID with negligible shielding, for which (28)-(29) reduce to φ Even in this simple case one can have FM coupling, for cos 2 While one thus needs a nonzero flux bias, tunability of I x b suffices for changing the sign of the coupling. One notes a few features, and it should be investigated whether these carry over to more generic devices described by our equations: (a) to achieve FM coupling, one needs nonzero flux and current bias; (b) the denominator in (31) is always positive; (c) cos φ (0) b1 and cos φ (0) b2 never become negative simultaneously. Need for current bias: for a symmetric SQUID this is easy to prove, since for I x b = 0 one then has φ b2 , and for any φ x b there is a stationary state with |φ (0) b1 | ≤ π/2. This cannot be generalized further, since we have already seen that the (asymmetric) dc-SQUID contains the rf-SQUID as a special case, and the latter can mediate FM coupling without any current bias. Need for flux bias: again easily proved for symmetric devices with arbitrary shielding, which have φ b2 etc. for φ x b = 0; in general, however, (28)-(29) show that an inductance imbalance between the two SQUID arms can play the same role as a nonzero φ x b . Finally, note that the rhs in (30) features the second-derivative matrix for the potential of the free b-SQUID (characterizing the device's susceptibility). For a stable minimum, this matrix should be positive-definite, from which (b) and (c) immediately follow in general. (As remarked above, the second derivatives keep appearing because of the M-expansion. For finite M, one has second differences instead [5].) One of the possible uses for a dc-SQUID coupler could be to moonlight as a readout device through switching-current measurement, so let us treat that in the present notation. Rather than think of the maximum bias for which (28)-(29) have solutions, it is convenient to fix φ x b and consider I x b as a function of the junction phases, which are still related among themselves by one nontrivial constraint. Critical bias now corresponds to dI x b /dφ = 0; working things out readily yields as the relevant condition. 3 To calculate switching-current sensitivity, we now vary φ x b , and see which variation in I x b conserves the above relation. One finds using (32), one again finds that this is more symmetric in the SQUID arms than it looks. So the sensitivity vanishes iff, at critical bias, both junctions become critical individually. This happens for As long as the "critical flux" productL bj I c,bj is different for the two arms, one has a finite sensitivity at φ x b = 0 and in fact for any φ x b , since one can ramp I x b with either sign. For the sake of completeness we also give the zero-sensitivity points which correspond to minima in the switching current. Finally, the above analysis only shows that the potential part of (24)-(25) enables tunable coupling in the classical case. That is, there could be devices outside its regime of validity which are nonetheless suitable couplers, as long as the transformer remains passive, i.e., confined to its lowest energy state/band [5]. For a quantum analysis, the interacting, biased Hamiltonian (24)-(25) can be used as a starting point. The method of M-expansion should again reduce the problem to one for the uncoupled, biased, dc-SQUID. Below, however, the quantum case is taken up for rf-SQUIDs only.

III. QUANTUM ANALYSIS
Given that quantum tunneling is crucial for the qubits a and c, a consistent quantum treatment of the whole system is desirable. In the quantum case, the M-expansion works as above. The other crucial step is an adiabatic approximation for the b-coupler (cf. [5]), replacing operators with their ground-state expectation. This very language is only meaningful if the rf-SQUID b is a well-defined separate entity, i.e., for weak coupling. Specifically, the coupling energies must be much smaller than the first excitation energy of the b-device, but this does not restrict their size compared to the level splittings of qubits a and c.
Thus, let us collect all terms in the Hamiltonian H = j=a,b,c Q 2 j /2C j − U [with U as in (1)] involving the coupler. To the relevant order in M one has with The adiabatic step amounts to the replacement and for small coupling, one can further write Here, is the ground-state expectation of the loop current; the susceptibility is defined as which gives the bound χ All that remains is to substitute (40) and (41) back into the full H; after some cancellation, one finds This does not involve assumptions on the states of the a-and c-qubits, or on their flux biases. Note that (46) involves current operators for the a-and c-devices [10], in contrast to, e.g., (42) and (44). Equations (44) and (45) imply that the qubits' effective parameters depend subtly on the coupler's working point, presenting a possible experimental challenge. As in the classical case, the perturbative step has reduced the problem to that of the uncoupled b-SQUID, and we henceforth omit its index as well as the zero superscript indicating free devices. Evaluating (42) in the potential minimum, one verifies that (46) reduces to (9) in the classical limit. In the quantum case, numerical differentiation in (42) presents no problems, but it will be instructive, and useful for Sec. IV, to express χ in terms of quantities at a single bias point. To this end, (42) is evaluated in terms of ξ 0 (φ) ≡ ∂ φ x ψ 0 (φ), where normalization implies that ψ 0 |ξ 0 = 0. In its turn, ξ 0 is solved from the φ x -differentiated Schrödinger equation in which ψ 0 figures as an inhomogeneous term, and where Using a variation-of-constant solution to (47), one then derives 4 This already expresses considerable simplification, in that two of the integrals in a triple integral were written as the square of will not be realized numerically for the former (latter) representation of k. Hence, we have to use both. Putting . Now, χ can be found by a single integration of the ODE systems A different derivation uses that, for stationary states, the shielding current equals the Josephson current (readily proven formally): It is correct to evaluate (50) with the aid of (47), but this yields unwieldy asymmetric expressions. Instead, we observe that a change in flux bias is a uniform translation if there is no Josephson potential, i.e., we set ξ 0 ≡ η 0 − ψ ′ 0 with ψ 0 |η 0 = 0 and using the matching expression E ′ 0 = E sin φ 0 . Combination leads to The translation of (52) to ODE form is wholly analogous to (49) and will not be repeated. For small β, (52) is clearly preferable to (48), since the double integral yields only a small correction to the first term. Compare (50), readily yielding an estimate of χ for weakly shielded SQUIDs, to (42), where in this case one has to evaluate the difference of two nearly equal terms. Conversely, for large β, φ 0 ≈ 2πn is approximately constant over most of the bias range, giving a clear advantage to (42) and (48), which capture this behaviour in their first terms. The equivalence of the four expressions for χ has been thoroughly tested numerically.
Results are shown in Fig. 4. It is seen that the finite spread of the ground-state wave function smooths out the classical susceptibility, in particular near φ x = π where the FM response is maximal. The finite quantum tunneling amplitude also means that there is no longer a sharp transition at β = 1, although the small tunnel splitting for larger β is likely to violate the adiabaticity requirement.

A. Analytics
While the above was derived with monostable coupling devices in mind, it should be general, hence equally applicable to rf-SQUID qubits close to the anticrossing, where it can be compared to the predictions of the two-state model with δΦ x ≡ Φ x − 1 2 Φ 0 . The bias ǫ is defined as half the difference between the "local well energies", which for definiteness are taken as the ground-state energies with an artificial nodal condition imposed at the potential maximum. If β − 1 is taken fixed, then the twostate model should be valid asymptotically in g ≡ 2CE/e 2 . Note that it incorporates several quantum effects short of tunneling: at finite bias the curvatures of the two wells change in opposite directions, affecting the zero-point energies (note 15 in Ref. [11]), and the anharmonicity of the wells is also accounted for. The "tunneling amplitude" ∆ is primarily defined through (53), without reliance on WKB estimates. Bearing in mind that ∆(δΦ x ) [ǫ(δΦ x )] is even [odd], one finds χ qb (0) = E qb 0 (0) ′′ = −∆ ′′ (0) − I 2 p /∆ for the two-state approximation to χ, with ∆ ≡ ∆(0) and I p ≡ −ǫ ′ (0).
These predictions make it clear that χ(0) is exponentially large in g, so that in the evaluation of (48) and (52), their first terms can be safely ignored. Further, φ 0 = π and sin φ 0 = 0 by symmetry for φ x = π, which will be assumed henceforth. In either formula, the increase in χ can only come from the factor ψ −2 0 , which is large for |φ| → ∞ and near the barrier. However, in both formulas and both for φ → ∞ and φ → −∞, this is overcome by {· · · } 2 decaying as ψ 4 0 . Therefore, the only relevant contribution comes from the barrier region, where to within exponentially small corrections the factors {· · · } 2 are constant, and readily related to the loop current. Thus, in the regime of interest, both (48) and (52) evaluate to where b means an integral over the barrier region, which for definiteness can be taken as everything between the two well minima (so as not to rely on "classical turning points" this defines h only up to a constant, since f was defined only up to a, possibly E-dependent, normalization. We now fix π −∞ dφ f 2 (φ) = 1 2 , so that h ′ (φ) ≈ −C/4e 2 f 2 (φ) throughout the barrier region. Further, with φ min the left potential minimum, h(φ min , E) is of order one, hence negligible, and one finds h(π, , and substitution in (54) yields Comparing with what was found below (53) for χ qb (0) in the two-state approximation, it leads one to conclude that actually ∆ ′′ (0) = 0. More precisely, comparing the orders of magnitude, it indicates that ∆ does not have relative variations of order one in the interval |ǫ| ∼ ∆, which is all one cares about. This flatness of ∆ in the anticrossing region is often assumed without much discussion, let alone proof; even the comparatively detailed WKB treatment in Ref. [12] is complete only for highly excited states in the local wells. Yet, the standard WKB expression for the tunneling amplitude is very sensitive to the barrier action integral, suggesting that there could be an equally strong sensitivity to bias variations. All in all, the above paints a positive picture of the standard two-level approximation. Its regime of validity should follow from |ǫ(δΦ x )| ≪ E 2 − E 1 , where the rhs is of the order of the plasma frequency. Hence, the range of allowed δΦ x is actually algebraically decreasing in g so that asymptotically, variations in loop current etc. become negligible within this range. At the same time, the range of permissible ǫ/∆ increases exponentially, so that wide sweeps in ǫ, as e.g. studied in the context of Landau-Zener transitions [13], can be performed without leaving this range.

B. Numerics
If one numerically finds the two lowest eigenfunctions of a near-degenerate bistable rf-SQUID, one can uniquely determine the parameters of the corresponding two-level lowenergy Hamiltonian 5 H qb = −ǫσ z − ∆σ x .
(57) 5 We have used similar wordings to describe the procedure in Refs. [3,11], but the concepts are different.
There, the impedance response as a function of flux bias was used to fit a single set of qubit-Hamiltonian coefficients. Here, we have numerical access to the wave functions, and can therefore determine these coefficients independently at each point in parameter space.
Results are summarized in Table I. With increasing β, one sees the same trends predicted at the end of Sec. IV A for increasing g: the range of validity of (57) with constant ∆ increases, in terms of not δΦ x but ǫ/∆. Remaining deviations ∆(δΦ x ) − ∆(0) ∼ (δΦ x ) 2 are too small to contradict the analytical conclusions. The breakdown of (57) occurs when the wave function in the depopulated well is so small that contributions near φ b can no longer be neglected, precluding a coarse-grained description. Even then, the linearity of ǫ(δΦ x ) is still excellent, as shown by the last entries in Table I for each β, and the breakdown is of little practical consequence as ∆ is too small compared to |ǫ| to have an appreciable effect on the spectrum of (57).
The above situation for a single qubit is uncharacteristically fortunate, in that the requirement of unitarity on U has circumvented coarse-graining the excited state ψ 1 . Already for two qubits, there is no natural counterpart to this shortcut, and the 4×4 matrix U having coarse-grained eigenstates for its columns is not exactly unitary, introducing a small uncertainty into the above procedure. The root cause of these subtleties is the problematic nature of the ubiquitous term "flux basis," which strictly speaking should consist of δ-functions in the coordinate representation. Pending further research, a completely satisfactory resolution may be as follows: calculate the matrix elements of φ in the low-energy subspace. This Hermitian matrix has an eigenbasis, and the eigenvectors can be taken as columns of a transition matrix V . 7 Recognizing that V consists of φ-eigenvectors in the energy basis whereas U consists of E-eigenvectors in the flux basis, one should thus consider H qb ≡ V † H E V as an alternative to (59).

V. CONCLUSION
A different sign-tunable flux transformer, relying on the gradiometric nature of the employed qubit, was previously considered in [14]. The transformer itself is also gradiometric, and tunability is achieved by incorporating compound junctions with variable couplings in either arm. In another type of gradiometric flux transformer, each arm of the gradiometer couples to one adjacent device, and the tunable element is a single compound junction in the central leg [15]. This does not give sign-tunability, though: compare (4) in [15] to (10) above. For the present type of mediated coupling, the prediction of sign-tunability is extremely generic: for any SQUID-type device, the I (0) (Φ x ) curve is smooth and periodic, and therefore will have pieces with both positive and negative slope in any Φ 0 -interval. Comparing the two implementations, the rf-SQUID (Sec. II A) does not have any galvanic coupling to external circuitry, which should limit decoherence. On the other hand, it was already mentioned that a dc-SQUID coupler (Sec. II B) [6] could be useful for readout. Either seems promising, and we look forward to experimental results on tunable coupling, which may well be the next big step forward in superconducting qubit technology.
The analysis in Sec. IV is a prelude to the case of multiple interacting qubits, where the precise form of the low-energy Hamiltonian H qb is still open to discussion. Namely, there could be multi-qubit tunneling, say, along the diagonals in a two-qubit potential, which would lead to σ x -σ x and σ y -σ y interactions ( [10], note in proof). There is little doubt that such terms would be small, but this does not automatically imply that their effect is negligible. Namely, they could, say, cause a gap at an |↑↓ ↔ |↓↑ anticrossing (for AF coupled qubits) in first order, while conventional tunneling terms ∼ σ x contribute to this gap only in second order. Further, there seems to have been little investigation of the possibility that (single-qubit) tunneling in qubit a might depend on the state of qubit b and vice versa, which would lead to σ x a σ z b etc. terms. To the extent that, for a-tunneling, the b-qubit can be regarded as an external flux bias, the flatness of ∆ observed in Sec. IV indicates that such terms should be very small. Finally, for more than two qubits, one can anticipate three-qubit etc. interactions-particularly if, for stronger coupling, one qubit's persistent current can depend on the state of the others. In any case, rotation methods such as in Eq. (59) should enable one to settle all these questions conclusively.