Qubit lattice coherence induced by electromagnetic pulses in superconducting metamaterials

Quantum bits (qubits) are at the heart of quantum information processing schemes. Currently, solid-state qubits, and in particular the superconducting ones, seem to satisfy the requirements for being the building blocks of viable quantum computers, since they exhibit relatively long coherence times, extremely low dissipation, and scalability. The possibility of achieving quantum coherence in macroscopic circuits comprising Josephson junctions, envisioned by Legett in the 1980’s, was demonstrated for the first time in a charge qubit; since then, the exploitation of macroscopic quantum effects in low-capacitance Josephson junction circuits allowed for the realization of several kinds of superconducting qubits. Furthermore, coupling between qubits has been successfully achieved that was followed by the construction of multiple-qubit logic gates and the implementation of several algorithms. Here it is demonstrated that induced qubit lattice coherence as well as two remarkable quantum coherent optical phenomena, i.e., self-induced transparency and Dicke-type superradiance, may occur during light-pulse propagation in quantum metamaterials comprising superconducting charge qubits. The generated qubit lattice pulse forms a compound ”quantum breather” that propagates in synchrony with the electromagnetic pulse. The experimental confirmation of such effects in superconducting quantum metamaterials may open a new pathway to potentially powerful quantum computing.


I. HAMILTONIAN FUNCTION AND QUANTIZATION OF THE QUBIT SUBSYSTEM
Consider an infinite number of superconducting charge qubits in a transmission line (TL) that consists of two superconducting plates separated by distance d, while the center-to-center distance between qubits, , is of the same order of magnitude. The charge qubit is of the form of a mesoscopic superconducting island which is connected to each electrode of the TL with a Josephson junction (JJ). Assume that an electromagnetic (EM) wave corresponding to a vector potential A = A z (x, t)ẑ propagates along the superconducting TL, in a direction parallel to the superconducting electrodes and perpendicular to the direction of the EM wave propagation. A minimalistic description of that system constists of writing the Hamiltonian of the compound qubit array -EM field system with the geometry shown in Fig.  1 of the paper, as H = n φ 2 n − 2 cos α n cos ϕ n +α 2 n + β 2 (α n+1 − α n ) 2 , where ϕ n is the superconducting phase on n−th island, β 2 = (8πdE J ) −1 (Φ 0 /(2π)) 2 , with d being the separation between the superconducting electrodes of the TL, and the overdots denote derivation with respect to the temporal variable t. In what follows, a basic assumption is that the wavelength λ of the EM field is much larger than the other length scales determined by the size of the unit cell of the qubit metamaterial (λ >> , d), so that the vector potential component A z (x, t) can be regarded to be approximatelly constant within a unit cell, A z (x, t) A z,n (t). Then, the discretized and normalized EM vector potential at the n−th unit cell which appears in Eq. (1) reads a n (t) = (2πd/Φ 0 )A n (t). The Hamiltonian function Eq. (1) is given in units of the Josephson energy E J = (Φ 0 I c )/(2πC), where I c and C is the critical current and capacitance, respectively, of the JJs, and Φ 0 = h/(2e) is the flux quantum, with h and e being the Planck's constant and the electron charge, respectively. That Hamiltonian can be written in a more transparent form by adding and subtracting 2 cos φ n and subsequently rearranging to get where the qubit subsystem energy H QU B , the EM field energy H EM F , and their interaction energy H int , take respectively the form In the following, the EM field is treated classically, while the qubits are regarded as two-level systems. The latter approximation is particularly well suited in the case of resonant qubit -EM field interaction adopted here. The quantization of the qubit subsystem can be formally performed by replacing the classical variables ϕ n anḋ ϕ n by the quantum operatorsφ n andφ n → −i ∂ ∂ϕn , respectively, in the Hamiltonian H QU B since we are dealing with a large number of Cooper pairs. The exact energy spectrum E p (n) and the corresponding wavefunctions Ξ p (n) of the nth qubit may then be obtained by mapping the Schrödinger equation with the single-qubit Hamiltonian H sq =φ 2 n − 2 cos ϕ n , onto the Mathieu equation The second quantization of the qubit subsystem proceeds by rewritting the single-qubit Hamiltonian as whereΨ † andΨ are field operators. Note that the subscript "n" in Eq. (5) has been dropped since the qubits are identical. Using the expansionΨ(ϕ) = p a p Ξ p (ϕ), where the operators a † p (a p ) create (annihilate) qubit excitations of energy E p , the Hamiltonian Eq. (5) is transformed into We hereafter restrict H sq to the Hilbert subspace of its two lowest levels, i.e., those with p = 0, 1, so that in second quantized form the Hamiltonian Eq.
(2) reads H = n p E p (n)a † n,p a n,p + p,p V p,p (n)a † n,p a n,p sin 2 α n 2 + where p, p = 0, 1 and the V p,p (n) ≡ V p ,p (n) that represent the matrix elements of the nth qubit -EM field interaction, are given by In the reduced state space, in which a single qubit can be either in the ground (p = 0) or in the excited (p = 1) state, the normalization condition p a † n,p a n,p = 1 holds for each qubit in the metamaterial. The reduced Hamiltonian Eq. (7) is also written in units of E J .

II. MAXWELL-BLOCH FORMULATION OF THE DYNAMIC EQUATIONS
In accordance with the adopted semiclassical approach, the time-dependent Schrödinger equation in whichH is the Hamiltonian from Eq. (7) in physical units, i.e.,H = HE J , is employed for the description of the qubit subsystem. Generally, the state of each qubit is a superposition of the form for which it can be easily shown that the coefficients Ψ n,p satisfy the following normalization conditions p |Ψ n,p (t)| 2 = 1, where in the second Eqs. (11) a finite N −qubit subsystem has been implied. Assuming that the pulse power is not very strong, substantial simplification of the dynamic equations for the qubit subsystem can be achieved using the approximation [1 − cos(α n )] (1/2)α 2 n in the interaction Hamiltonian H int . Substitution of |Ψ = |Ψ n from Eq. (10) into the Schrödinger equation (9), along with derivation of the classical Hamilton's equation for the normalized EM vector potential α n , yields where χ =hω j /E J . In Eqs. (12) and (13), the temporal variable is renormalized according to t → ω J t. Because of that, the dimensionless energy of the qubit excitations has to be redefined according to E P → p = E p /χ, which explains the presence of the dimensionless factor 1/χ in Eq. (12).
in which the variables R i (i = x, y, z) apply to each single qubit. Thus, from the point of view of a macroscopic system, |Ψ i | 2 represent the population densities i = N N |Ψ| 2 ≡ Ni N . By transforming Eqs. (12) and (13) according to Eq. (14) we geṫ , By taking the continuous limit of Eqs. (15) and (16) using the relations α n → α(x, t) and α n±1 ≈ α ± α x + 1 2 α xx , we obtainṘ where R x , R y , R z , and α are functions of the spatial variable x and normalized temporal variable t, while the overdots denote partial derivation with respect to the latter. It can be easily verified that the Bloch equations Eqs. (18) possess the dynamic invariant i R 2 i = 1, whose constant value is derived through the normalization conditions.

III. SLOWLY VARYING ENVELOPE APPROXIMATION (SVEA)
The Slowly Varying Envelope Approximation (SVEA) relies on the assumption that the envelop of a travelling pulse in a nonlinear medium varies slowly in both time and space compared with the period of the carrier wave, that makes possible the introduction of slow and fast variables. Assuming that in the following, we can approximate the EM vector potential function by where ψ(x, t) = kx − ωt + φ(x, t), with k and ω being the wavenumber and frequency of the carrier wave, respectively, that are connected through the dispersion relation that is obtained at a later stage. The functions ε(x, t) and φ(x, t) in Eq. (20) are the slowly varying envelop and phase, respectively. Using fast and slow variables, the x and y Bloch vector components, R x (n) and R y (n), can be expressed as a function of new, in-phase and out-of-phase Bloch vector components r x and r y as A. Derivation of the dynamic equations for the slowly varying envelop and phase of the electromagnetic vector potential From Eqs. (20) and (21), the second temporal and spatial derivative of α(x, t) can be approximated bÿ in which the rapidly varying terms of the formε, ε xx , φ 2 , φ xx ,φ, φ x ε x , etc., have been neglected. Substitution of Eqs. (22) and (21), into Eq. (19) gives 2(ωε + kβ 2 ε x ) sin ψ + [2(φω + kφ x ) − ω 2 + Ω 2 + β 2 k 2 ]ε cos ψ = −χ{Dr z + µ[r x cos(2ψ) + r y sin(2ψ)]}ε cos ψ. (23) Equating the coefficients of sin ψ and cos ψ in the earlier equation yields and In order to obtain the dispersion relation ω = ω(k), we require the vanishing of the expression in the curly brackets in Eq. (25), which yields for the wavevector of the EM radiation (i.e., the pulse) in the TL the expression Thus, the EM wave can propagate through the superconducting quantum metamaterial only when its frequency exceeds a critical one, ω c = Ω = (V 00 + V 11 )/2. Finally, Eqs. (24) and (25) are averaged in time over the period T = 2π/ω (i.e., the fast time-scale) of the phase ψ(x, t). Due to the simple time-dependence of ψ(x, t) that has been assumed earlier in the framework of SVEA, that averaging actually requires the calculation of integrals of the form This procedure, when it is applied to the two evolution equations Eqs. (24) and (25) gives the truncated equations for slow amplitude and phaseε where c = β 2 k/ω is a critical velocity.

B. Derivation of dynamic equations for the transformed Bloch vector components
Substituting Eqs. (20) and (21) into Eqs. (18) for the original Bloch vector components, we get In order to derive the dynamic equations for the new Bloch vector components r i , we first multiply Eqs. (28) and (29) by cos(2ψ) and sin(2ψ), respectively, and then subtract one equation from the other. Thus we geṫ r x + 2ψr y = −(∆ + 2Dε 2 cos 2 ψ)r y + 2µε 2 cos 2 ψ cos(2ψ)r z .
The truncated Bloch equations Eqs. (33), along with Eqs. (27) for the slowly varying amplitude and phase of the EM vector potential pulse, constitute a closed system with five unknown functions which describe the approximate dynamics of the superconducting quantum metamaterial that are determined in the next section.

IV. EXACT INTEGRATION OF THE TRUNCATED EQUATIONS
The combination of Eqs. (27) and the third equation of Eqs. (33) provides a relation between the slow amplitude and the phase of the EM vector potential pulse. The first of Eqs. (27) is multiplied by ε and written as while the derivative of the second of Eqs. (27) is taken with respect to time and thenṘ z =ṙ z is replaced from Eqs.
(33). Thus we get From Eqs. (35) and (36), and by taking into account the independence of the slow temporal and spatial variables, we get where the constant of integration can be set equal to zero. Using Eq. (37), the truncated Bloch equations Eqs. (33) can be written asṙ The latter can be written in a simpler form using the unitary transformation where Φ is the constant transformation angle that is going to be determined. The truncated Bloch equations for the r i , i = x, y, z, can be written in terms of the new Bloch vector components S i , using a procedure similar to that used in Subsection III.B to obtain Eqs. (33). Sustituting Eqs. (39) into Eqs. (38), we geṫ Multiplying Eqs. (40) and (42) by cos Φ and sin Φ, respectively, and then adding them together, we geṫ Similarly, multiplying Eqs. (40) and (42) by sin Φ and cos Φ, respectively, and then subtracting one equation from the other, we getṠ The transformation angle Φ can now be selected so that the resulting equations are simplified as much as possible. We thus define so that cos Φ = ±σ and sin Φ = ±σγ where σ = 1/ 1 + γ 2 . The choice of the sign is irrelevant and here we pick positive sign for both functions. Using the above value for the transformation angle and the definitions W = (4D) 2 + µ 2 and η = −δµ/W , Eqs. (43), (41), and (44) obtain their final forṁ In order to investigate the possibility for "coherent propagation" of an electromagnetic potential pulse, we consider resonance conditions. In that case, η = 0, and Eqs. (46) becomė Combining the second and third Eq. (47) and integrating, we obtain the "resonant" conservation law S 2 y +S 2 z = const.. Assuming that all the qubits are in the ground state at t = −∞, we have the initial conditions r x (t = −∞) = r y (t = −∞) = 0 and r z (t = −∞) = −1 which are transformed into S x (t = −∞) = −γσ, S y (t = −∞) = 0, and S z (t = −∞) = −σ through the transformation Eq. (39). Applying the initial conditions to the resonant conservation law, we get In the following, we seek solution of the form ε = ε(τ = t − x/v) and S i = S i (τ = t − x/v), with i = x, y, z. By changing the variables in the first of Eqs. (27), with r y being replaced by S y , we get after rearramgement Then, combining Eq. (29) with the third Eq. (47) and integrating, we get where the conditions ε(−∞) and S z (−∞) = −σ were used. The system of Eqs. (48)-(50) for ε, S y , and S z can be integrated exactly; the variables S y and S z can be eliminated in favour of ε to give ε τ = λε 2 √ a + bε 2 , in which the constants are defined as a = 2σ/κ ωW v c−v to simplify the notation. The equation for ε can be readily integrated where we have set ε 0 ≡ ε(τ = τ 0 ) = √ 2σκ to eliminate the boundary term resulting from the integral over ε. Solving Eq. (51) for ε, we finally get a Lorentzian-like slowly varying amplitude where

V. DETAILS OF THE NUMERICAL SIMULATIONS AND THE ROLE OF DECOHERENCE
The numerical simulations in order to confirm (or not) the existence of the quantum coherent effects and the quantum coherence induced in the medium, i.e., the qubit transmission line, which are predicted through analytical considerations, consist of integrating Eqs. 3a-3c and 4 of the paper (n = 1, N ) andα n + χ Ω 2 + µR x (n) + DR z (n) sin α n = β 2 δa n , where δα n = α n−1 − 2α n + α n+1 and the parameters D, Ω, p , µ, ∆, and χ are defined in the paper as functions of the interaction matrix elements V ij (i, j = 1, 2), the energy difference between the ground and the first excited level, E 1 − E 0 , and the discrete vector potential coupling parameter β. The earlier system of ordinary differential equations is integrated in time with a standard fourth order Runge-Kutta algorithm with constant time-step ∆t, typically 10 −3 . Such small time-steps, and even smaller sometimes, are required to conserve up to high accuracy the total and partial probabilities as the compound system of qubits and the electromagnetic vector potential evolve in time. In all the simulations, periodic boundary conditions are used. Due to the particular shape (Lorentzian) of the EM vector potential pulse and the population inversion pulse in the qubit subsystem, very large systems with N = 2 13 = 8192 and N = 2 14 = 16384 have been simulated to diminish as much as possible the effects from the ends (i.e., the interaction of the pulse tail with itself, as it could be in the case of periodic boundary conditions). In some cases, it is necessary to simulate even larger systems, with N = 50, 000. In the figures presented in the paper and here only part of these arrays is shown, for clarity. In order to observe two-photon self-induced transparent (TPSIT) pulses a n (t) and the induced quantum coherent population inversion pulses R z (n; t), the following initial coditions are implemented: for the former, the analytically obtained solution resulting for the given set of parameters, while for the latter all the qubits are set to their ground, E 0 , state. In terms of the Bloch vector variables R i , i = x, y, z, that initial condition is specified as: for any n = 1, ..., N . In that case, the above mentioned pulses exist for velocities less than the corresponding limiting velocity for TPSIT media, i.e., The last equality is valid only when the two-photon resonance condition ω = ∆/2 is satisfied. In Eq. (8), k and ω denote the wavenumber and frequency of the carrier wave of the electromagnetic vector potential; while the latter is determined by the two-photon resonance condition, the former may vary within an interval which is restricted by the qubit-parameter-dependent condition that guarantees the wavenumber k being real. The role of decoherence in TPSIT in superconducting quantum metamaterials like the one investigated here can be observed in Figures 1, 2, and 3, in which the amount of decoherence, quantified by the factor γ, takes the values 0, 0.01 and 0.1, respectively. The decoherence factor γ depends solely on the interaction matrix elements V ij ; it is given by i.e., it is equal to zero for V 00 = V 11 while it assumes non-zero values for V 11 > V 00 . The Figures 1-3 show snapshots of R z (n; t) and a n (t) at several instants from t = 0 to t = 168, which are separated by 14 time units, along with the corresponding analytical solutions. All these profiles but the first are shifted downwards to avoid overlapping, while only part of the array is shown for clarity. Note that time increases downwards. In Figure 1, for γ = 0, the a n (t) pulse is observed to excite a population inversion pulse R z (n; t) in the qubit subsystem ( Figure 1a). The amplitude of R z (n; t) gradually increases until it attains its maximum value close to unity, while at the same time it propagates to the right with velocity v . In the figure, that occurs for the first time at t 70 time units; subsequently it evolves in time while it keeps its amplitude almost constant for at least 56 time units. After that, its amplitude starts decreasing until it is completely smeared (not shown). During the time interval in which the amplitude of the R z (n; t) pulse is close to the predicted one (i.e., close to unity), the metamaterial is considered to be in the almost quantum coherent regime. Note that at about t = 84 a little bump starts to appear which grows to a little larger ("probability bump") which gets pinned at a particular site at n ∼ 300. Note that another such bump appears at n = 0 due to the initial "shock" of the qubit subsystem because of the sudden onset of the a n (t) pulse. A comparison of the numerical R z (n; t) profiles with the analytical ones reveals that the velocity of propagation v , which is the same as that of the numerical a n (t) pulse (Figure 1b), is slightly larger than the predicted one v (v > v). The corresponding profiles for very weak decoherence, in which the value of the decoherence factor γ = 0.01, are shown in Figure 2. For that weak decoherence, there are only slight differences in comparison with Figure 1 which actually cannot be observed in the scale of the figures.
In Figure 3, the decoherence factor has acquired a substantial value of γ = 0.1, so that its effects are now clearly observed in the population inversion pulse R z (n; t), while the vector potential pulse a n (t) does not seem to be affected. In Figure 3a, while the R z (n; t) is still excited by the vector potential pulse a n (t), it has a very low amplitude compared with that of the analytically predicted form. Indeed, that amount of decoherence only slightly changes the analytical propagating pulses in superconducting quantum metamaterials (SCQMMs) without decoherence. a, Snapshots of the population inversion pulse Rz(n; t), excited by the induced quantum coherence in the qubit subsystem by the electromagnetic vector potential pulse, in the absence of decoherence (γ = 0); the pulse propagates to the right (time increases downwards). b, Snapshots for the corresponding evolution of the electromagnetic vector potential pulse an(t), that exhibits significant broadening while its amplitude decreases as time passes by; the numerical and analytical pulses are shown in blue and red color, respectively. Parameter values: χ = 1/4.9, β = 6, V00 = V11 = 1 (γ = 0), V01 = V10 = 0.7, E1 − E0 = 3, and v/c = 0.7. R z (n; t) pulse and that change is hardly visible in the scale of the figure. Note that the speed of the R z (n; t) is the same as in the case without decoherence (which is also the same with that of the a n (t) pulse, in Figures 1-3); even the unwanted "probability bumps" appear at about the same locations with almost the same amplitude and shape independently of the amount of decoherence.
Next, the possibility for the existence of propagating two-photon superradiant pulses in SCQMMs for large values propagating pulses in superconducting quantum metamaterials (SCQMMs) in the presence of weak decoherence. a, Snapshots of the population inversion pulse Rz(n; t), excited by the induced quantum coherence in the qubit subsystem by the electromagnetic vector potential pulse, in the absence of decoherence (γ = 0.01); the pulse propagates to the right (time increases downwards). b, Snapshots for the corresponding evolution of the electromagnetic vector potential pulse an(t), that exhibits significant broadening while its amplitude decreases as time passes by; the numerical and analytical pulses are shown in blue and red color, respectively. Parameter values: χ = 1/4.9, β = 6, V00 = 0.998, V11 = 1.002, V01 = V10 = 0.7, E1 − E0 = 3, and v/c = 0.7.
of β is numerically explored (Figure 4). In order to obtain that figure, the qubit subsystem is initialized with all the qubits in their excited state, while the vector potential pulse assumes its analytically predicted form for the selected parameter set. A very large array with N = 50, 000 is simulated, and the initial position of the (center of) a n (t) pulse is at n = −18, 750. The subsequent evolution produces the snapshots presented in Figure 4, from t = 0 up to t = 300 FIG. 3: Numerical validation of the analytical expressions for two-photon self-induced transparent (TPSIT) propagating pulses in superconducting quantum metamaterials (SCQMMs) in the presence of substantial decoherence. a, Snapshots of the population inversion pulse Rz(n; t), excited by the induced quantum coherence in the qubit subsystem by the electromagnetic vector potential pulse, in the absence of decoherence (γ = 0.1); the pulse propagates to the right (time increases downwards). b, Snapshots for the corresponding evolution of the electromagnetic vector potential pulse an(t), that exhibits significant broadening while its amplitude decreases as time passes by; the numerical and analytical pulses are shown in blue and red color, respectively. Parameter values: χ = 1/4.9, β = 6, V00 = 0.98, V11 = 1.02, V01 = V10 = 0.7, E1 − E0 = 3, and v/c = 0.7. time units (time increases downwards). These snapshots are separated in time by 20 time units and they are displaced vertically to avoid overlapping. The population inversion R z (n; t) and the electromagnetic vector potential pulse a n (t) are shown in Figures 4a and 4b, respectively. In that simulation, the parameter β is much larger than that used in Figure 3b of the paper, i.e., here β = 30. The other parameters, except the pulse speed, here v = 2.5 > c, are the