Physical device and control

The studied system consists of an array of 15 twin flux qubits, each of them containing five Josephson junctions (see Fig. 1). The central Josephson junction (called α-junction) of each qubit is shared by two superconducting loops and embedded directly into the central electrode of an Al coplanar waveguide (CPW). This setup provides an extremely large coupling between qubits and electromagnetic waves propagating through the system. We measured the microwave transmission through the CPW as a function of frequency and an external dc magnetic field applied perpendicular to the plane of the substrate. Experiments were performed at low power of the applied microwave field to remain in the linear driving regime. The radiation power is estimated to be on the order of −125 dBm at the chip input port.

Fig. 1 A quantum metamaterial made of twin flux qubits. Superconducting quantum metamaterial consisting of an array of 15 twin qubits embedded in a coplanar wave guide. An SEM image of twin flux qubits (above) and a whole structure (below) are shown. Each qubit consists of two superconducting loops sharing one common central Josephson junction (α-junction) and four identical Josephson junctions located on the outer parts of the loops. The α-junction allows the magnetic flux to tunnel between the loops. The inset is a schematic of a single meta-atom—the twin flux qubit; the phases on nodes are shown Full size image

Twin flux qubit model

Precluding the following demonstration of the experimental data and its detailed theoretical analysis, we present the essence of our idea to implement a metamaterial based on twin flux qubit. The twin flux qubit circuit has four nodes. By grounding one of the nodes (setting its nodal phase to zero) and assigning nodal phase variables to the rest of them we can completely describe the state of the system. The Hamilton’s function, which is sum of the charge and Josephson energies, reads

$$\begin{array}{l}H(\phi ,{{q}}) = \frac{{(2e)^2}}{2}{{q}}\hat {{C}}^{ - 1}{{q}}^T + E_{\mathrm{J}}\left( {4 + \alpha - \alpha \,{\mathrm{cos}}\,\phi _0} \right.\\ \left. { - {\mathrm{cos}}\,\phi _1 - {\mathrm{cos}}\,\phi _2 - {\mathrm{cos}}\left( {\phi _2 - \phi _0 + \phi } \right) - {\mathrm{cos}}\left( {\phi _1 - \phi _0 - \phi } \right)} \right)\end{array}$$ (1)

where \(\phi = \left( {\phi _0,\phi _1,\phi _2} \right)\) are the nodal phases and q = (q 0 , q 1 , q 2 ) are the corresponding canonically conjugate charges (see inset Fig. 1). Here, Φ is the external flux applied to a single loop of the twin qubit. Since the Josephson energy is periodic in Φ, we analyze the phase distribution across the twin qubit structure at applied fluxes Φ from 0 to Φ 0 per single loop of the twin-qubit, where Φ 0 = h/2e is the magnetic flux quantum. When Φ is close to 0, the phase across the central junction ϕ 0 is zero. Simple calculation demonstrates that at Φ close to Φ 0 /2 the minimum energy corresponds to ϕ 0 = π. This result follows from phase quantization of the superconducting order parameter around the loop. Therefore, one can expect a sharp change of the transmission coefficient \(\left| {S_{21}} \right|^2\) at the critical value of Φ cr (between 0 and Φ 0 /2), and a jump back at the external flux \({{\Phi }} \simeq {{\Phi }}_0 - {{\Phi }}_{{{{\rm cr}}}}\). The value of Φ cr depends on circuits parameters.

Replacing the nodal phase and charge variables by the respective operators in Hamilton’s function (see Eq. (1)) yields the Hamiltonian operator. By making use of numeric diagonalisation of the Hamiltonian the eigenvalues and the eigenstates are found (see Supplementary Note 1). The magnetic field dependence of the ground state energy of a twin flux qubit is presented in Fig. 2a. For the chosen experimental parameters of the twin flux qubit, in the vicinity of zero-applied flux, the ground-excited state transition energy hf 01 is quite large (see Fig. 2b). However, if the external flux is close to the value of Φ 0 /2 per each loop, the energy hf 01 is about 10 GHz and has an extremum (“sweet spot”). At this point, the twin-qubit transition energy reaches a local maximum rather than minimum characteristic for the standard single-loop three-junction flux qubit20. Moreover, the flux dependence of f 01 is much weaker than that of the usual flux qubit. This makes the twin-qubit more immune to the well known flux qubit problem—strong sensitivity of the transition energy to the external magnetic field fluctuations leading to unavoidable dephasing.

Fig. 2 Magnetic field dependence of the transition frequencies of the twin qubit. The energy of the ground state (a) and the transition energy hf 01 of the twin qubit calculated from the Hamiltonian of Eq. (1) (b). The parameters α = 0.72 and C = 5.2 fF and the Josephson energy is E J = 50 GHz. These dependencies are Φ 0 periodic and symmetric with respect to Φ/Φ 0 = 0.5. The minimum point of the (b) plot corresponds to the transition of the central junction phase φ 0 from zero to π Full size image

Measurements of the Transmission Coefficient

The measured dependence of the transmission coefficient \(\left| {S_{21}} \right|\) on both the microwave frequency f and the externally applied magnetic flux is shown in Fig. 3a. First, as it is expected, one notices that \(\left| {S_{21}} \right|\) is symmetric around zero field and periodic with respect to the external magnetic flux Φ (outside the range of Fig. 3) with period of Φ 0 per single qubit loop. Second, a spectacular result of the measurement is the presence of two regions in the magnetic flux with drastically different transmission coefficients \(\left| {S_{21}} \right|\), see the orange and blue areas in Fig. 3a. The transition between these two regimes is abrupt and occurs at magnetic flux values Φ cr ≃ ±0.3Φ 0 , see Fig. 3c. In the range 0.3Φ 0 ≤ Φ ≤ 0.7Φ 0 , the transmission amplitude \(\left| {S_{21}} \right|\) changes weakly (from −10 to −5 dB) in frequency band 6–14 GHz outside the resonance (the red feature). The intriguing field-dependent resonant features are characterized by a large enhancement of the transmission amplitude (appearing in Fig. 3a in red) over a relatively narrow frequency range between 11 and 14 GHz. The two symmetric resonant features reach their maximum frequency near Φ ~ ±Φ 0 /2. The induced resonant transmission frequency is smoothly decreasing on either side from the maximum and disappears at around 11 GHz.

Fig. 3 Transmission of microwaves through the quantum metamaterial in different regimes. a The measured dependence of the amplitude of transmission coefficient t (normalized to the value at zero field) on applied dc magnetic field (proportional to the bias current in the coil, lower axis) and frequency f. The upper horizontal axis translates the field in magnetic flux Φ per qubit single loop. The transmission t displays the sharp changes under variation of the magnetic flux Φ. One can see two different ranges of microwave propagation, nearly flat transmission around zero field and sharp resonant enhancement of the transmission near 11–14 GHz at magnetic flux Φ ~ ±Φ 0 /2. b A cross-cut of a at the fixed frequency of 13 GHz. The sharp peaks correspond to coherent tunneling between quantum states in the twin qubits (see text). c A cross-cut of a at the fixed frequency of 10 GHz. The sharp jumps correspond to a transition between zero and π phase on the central junction of the twin qubit (see text). Red curve is a fit to the theoretically predicted dependence Eq. (12) Full size image

In order to understand qualitatively the observed modulation of the microwave transparency, let us consider the dependence of the qubit state on the applied magnetic field. At low temperatures, the metamaterial is in its ground state. Therefore, the phases across qubit junctions can be reconstructed as functions of magnetic fluxes applied through qubit loops. Since the external magnetic flux is applied symmetrically to the twin qubits structure, induced current across central junction is zero, and its phase, consequently, is either zero or π. Along with increasing the applied magnetic flux, a circulating current through external junctions is induced, providing a monotonic increase of the Josephson phases of the outer Josephson junctions and, in turn, the effective qubit inductance. This causes a smooth decrease of the transparency, see Fig. 3c.

As soon as the external magnetic flux reaches the value of Φ cr , it is energetically favorable to increase the energy of the central junction and to lower the energies of the outer junctions. Quantitatively this change is dictated by the generalized Bohr–Sommerfeld quantization rule ensuring a single-valued wave-function of the condensate across qubits loops. As a result at \(\left| {{\Phi }} \right| = {{\Phi }}_{{\mathrm{cr}}}\) the phase shift of the central junction ϕ 0 from zero to π is accompanied by the phase shift of the outer junctions by π/2, see Supplementary Fig. 4, and the transmission amplitude sharply drops, see Fig. 3c. We call the ground state of the twin qubit in this magnetic field range the π-state as opposed to the zero-state at \(\left| {{\Phi }} \right| < {{\Phi }}_{{\mathrm{cr}}}\). In the π-state, the currents in the outer branches of the twin flux qubit also change according to the new values of the phases. In particular, at the “sweet” spot Φ = Φ 0 /2 phases of all five Josephson junctions of the twin qubit are either zero or π which implies zero currents flowing through all junctions, causing the absence of current flowing in all branches of the twin flux qubit.

This nonresonant effect is observed for the entire frequency range used in our experiments and caused by the field-dependent phase distribution in the twin-qubit metamaterial in its ground state. The experimentally obtained value of Φ cr is consistent with the calculated dependence on α, see Supplementary Note 2. Important to note that widely discussed Electromagnetically Induced Transparency (EIT) has a different nature than the phenomenon described here. “Conventional” EIT is caused by the quantum interference between different scattering channels. It immediately demands the presence of many (at least three) energy levels for joint photon–qubit systems21,22,23,24,25,26. At contrary in the case considered here the classical electromagnetic field interacts with a regular array of qubits, and EIT occurs due to changing of the ground state of a twin flux qubits arrangement.

In the π-state region of external magnetic flux, the transition energy between the ground and first excited states of the twin qubit hf 01 can be determined experimentally. The sharp peaks in Fig. 3b and inverted red parabolas between 11 and 14 GHz in Fig. 3a correspond to resonant enhancement of transmission amplitude resulting from the \(\left| 0 \right\rangle \to \left| 1 \right\rangle\) transition. These experimental data in Fig. 3b are fitted by the transition frequency f 01 obtained from diagonalisation of the Hamiltonian (1) and excellent agreement is observed, also see Supplementary Fig. 2.

Analysis of the transmission coefficient

In spite of the quantum nature of the twin qubit and the strong coupling between qubits and applied microwaves, here we first describe the electromagnetic field as a classical wave. Below we quantitatively analyze the wave propagation through the transmission line which is strongly inductively coupled to a periodic array of quantum oscillators. The transmission line is characterized by coordinate and time dependent charge distribution Q(x, t). Here, x is the coordinate along the transmission line. In the presence of a periodic array of twin qubits, which can be considered as superconducting SQUID-like oscillators, the electromagnetic wave propagation is determined by the equation as

$$c_0^2\frac{{\partial ^2Q}}{{\partial x^2}} - \frac{{\partial ^2Q}}{{\partial t^2}} - \gamma _{{\mathrm{tl}}}\frac{{\partial Q}}{{\partial t}} + \mathop {\sum}\limits_n \beta \left\{ {Q\left( {x_n,t} \right)} \right\}\delta \left( {x - x_n} \right) = 0$$ (2)

where the parameter β{Q(x n , t)} is the effective scattering strength of a single twin qubit located at x = x n . The parameter γ tl describes the intrinsic dissipation of the transmission line (γ tl ≪ 1), and c 0 is the velocity of electromagnetic waves in the transmission line.

In the particular case of a twin flux qubit (see Fig. 1) the parameter β is obtained from the quantum-mechanical average of the time-dependent Josephson phase of the α-junction φ 0 (t), i.e.,

$$\beta = \frac{\hbar }{{2eL}}\left\langle {\dot \phi _0(t)} \right\rangle$$ (3)

where L is the transmission line inductance per unit length. The quantum dynamics of the Josephson phase ϕ 0 is determined by the stationary Hamiltonian (1) and the time-dependent perturbation Hamiltonian, \(H_{{\mathrm{per}}} = \left( {E_{\mathrm{J}}{\mathrm{/}}I_{\mathrm{c}}} \right)\dot Q(x,t)\phi _0\).

Let us consider an electromagnetic wave of a given frequency ω, Q(x,t) = Q(x)eiωt. In this case, Eq. (2) can be re-written as

$$c_0^2\frac{{\mathrm{d}^2Q(x)}}{{\mathrm{d}x^2}} + \left( {\omega ^2 - \gamma _{{\mathrm{tl}}}i\omega } \right)Q(x) + \mathop {\sum}\limits_n \beta \left\{ {Q\left( {x_n} \right)} \right\}\delta \left( {x - x_n} \right) = 0$$ (4)

Following the method elaborated for the solution of the Schrödinger equation with the Kronig–Penney potential27, we can present the charge distribution Q(x) in the following form:

$$Q(x) = - \frac{i}{{2k}}\mathop {\sum}\limits_n \beta \left\{ {Q_n} \right\}{\mathrm{exp}}\left( {ik\left| {x - x_n} \right|} \right)$$ (5)

where the wave vector \(k = \sqrt {\omega ^2 + i\gamma _{{\mathrm{tl}}}\omega } {\mathrm{/}}c_0\) and Q n = Q(x n ) is the amplitude of propagating charge distribution at point x n . By making use of the properties of the δ-function, we obtain a set of discrete equations for Q n :

$$Q_{n + 1} + Q_{n - 1} - 2\left[ {Q_n\,{\mathrm{cos}}\,ka - \frac{{\beta \left\{ {Q_n} \right\}}}{{2kc_0^2}}{\mathrm{sin}}(ka)} \right] = 0$$ (6)

Here, a is the effective distance between adjacent point-type scatterers of the array. Taking into account that, in the microwave range, \(ka \ll 1\), we substitute the discrete values of Q n with a smooth function \(\tilde Q(x)\) satisfying the equation

$$\frac{{{\rm d}^2\tilde Q(x)}}{{\mathrm{d}x^2}} + k^2\tilde Q(x) - \frac{{\beta \left\{ {\tilde Q(x)} \right\}}}{{ac_0^2}} = 0$$ (7)

For low amplitude of the waves, we obtain an explicit expression for the effective scattering strength β:

$$\beta = - \frac{{\hbar \omega ^2}}{{2eLI_{\mathrm{c}}}}K(\omega )Q,K(\omega ) = \frac{{\omega _{\mathrm{p}}^2}}{{{{\Omega }}^2 - \omega ^2 + i\gamma _{\mathrm{J}}}}$$ (8)

where Ω = ΔE/ħ is the characteristic frequency of the oscillator that can be used to approximate the lowest-lying energy levels of the twin flux qubit. For small external flux values \(\left| {{\Phi }} \right| < {{\Phi }}_{{\mathrm{cr}}}\), where Φ cr = Φ 0 /π arcsin α, the 0-state is the unique stable state of the twin qubit. For the qubits with α value used in this paper, Φ cr ≈ 0.3Φ 0 . The characteristic frequency Ω = Ω + in the range \(\left| {{\Phi }} \right| < {{\Phi }}_{{\mathrm{cr}}}\) is relatively large (see Fig. 2b).

$${{\Omega }}_ + ^2 = \omega _{{p}}^2\frac{{{\mathrm{cos}}\left[ {\pi {{\Phi /\Phi }}_0} \right] + \alpha }}{{1 + \alpha }}$$ (9)

Notice here, that the frequency Ω + smoothly decreases with the applied magnetic flux Φ and results in a parabolic non-resonant decrease of \(\left| {S_{21}} \right|\) as observed in experiments, see the central part of the \(\left| {S_{21}} \right|\) curve in Fig. 3b.

The characteristic frequency of small-amplitude quantum oscillations in the π-state for \({{\Phi }}_{{\mathrm{cr}}} < \left| \Phi \right| < {{\Phi }}_0 - {{\Phi }}_{{{{\rm cr}}}}\) can be written explicitly as (see the central part of Fig. 2b).

$${{\Omega }}_ - ^2 = \omega _{{p}}^2{\kern 1pt} \frac{{{\mathrm{sin}}\left[ {\pi {{\Phi /\Phi }}_0} \right] - \alpha }}{{1 + \alpha }}$$ (10)

Since the characteristic frequency \({{\Omega }}_ - \ll {{\Omega }}_ +\) one expects strong changes of the magnitude of the transmission coefficient \(\left| {S_{21}} \right|\) at magnetic flux values Φ = Φ cr and Φ = Φ 0 − Φ cr .

Substituting (8) in (7) we obtain an equation for the propagation of electromagnetic waves in this quantum metamaterial:

$$\frac{{{\rm d}^2\tilde Q(x)}}{{\mathrm{d}x^2}} + k^2[1 + AK(\omega )]\tilde Q(x) = 0$$ (11)

where the parameter A ≈ ħ/(2eLaI c ) depends on the properties of both the transmission line and the Josephson junctions. In our case \(A = 400 \gg 1\) and in this limit we can obtain the transmission coefficient \(S_{21} = \left| {Q(l){\mathrm{/}}Q(0)} \right|\) (l is the length of a transmission line) in the following form: