Chiral SQUID-metamaterial waveguide for circuit-QED

Superconducting metamaterials, which are designed and fabricated with structured fundamental circuit elements, have motivated recent developments of exploring unconventional quantum phenomena in circuit quantum electrodynamics (circuit-QED). We propose a method to engineer 1D Josephson metamaterial as a chiral waveguide by considering a programmed spatiotemporal modulation on its effective impedance. The modulation currents are in the form of traveling waves which phase velocities are much slower than the propagation speed of microwave photons. Due to the Brillouin-scattering process, non-trivial spectrum regimes where photons can propagate unidirectionally emerge. Considering superconducting qubits coupling with this metamaterial waveguide, we analyze both Markovian and non-Markovian quantum dynamics, and find that superconducting qubits can dissipate photons unidirectionally. Moreover, we show that our proposal can be extended a cascaded quantum network with multiple nodes, where chiral photon transport between remote qubits can be realized. Our work might open the possibilities to exploit SQUID metamaterials for realizing unidirectional photon transport in circuit-QED platforms.

In circuit-QED the boson modes of microwave photons are routinely supported by conventional circuit elements such as LC resonators, transmission-line resonators and 1D open coplanar-waveguides [29][30][31][32][33]. Recently exploring intriguing QED phenomena with metamaterial structures in SQCs have attracted a lot of interests [34][35][36][37][38][39][40]. Compared with standard circuit-QED elements, superconducting metamaterials can be designed and fabricated in various kinds of spatial structures, which might have exotic band spectrum and nontrivial vacuum modes [41,42]. For example, by engineering the hopping rates among microwave resonators or SQC qubits, the photonic metamaterials analogue to topological lattices and strongly correlated matters were realized [43,44]. When the inductors and capacitors in the usual discrete model of a 1D transmission line are interchanged, the * wangxin.phy@xjtu.edu.cn left-handed metamaterials can also be fabricated, which might work as a waveguide with an effective negative index of refraction [45][46][47][48][49].
For example, many-body quantum optics in ultra-strong coupling regimes was successfully demonstrated in a Josephson chain platform by exploiting its high characteristic impedance [59]. Moreover, given that each node junction is replaced with a superconducting quantum interference device (SQUID), the SQUID chain is possible to be engineered as a tunable high-impedance ohmic reservoir by applying position-dependent magnetic flux [60]. When the SQUID chain are biased with a spacetime varying flux, analogue Hawking radiation can be reproduced [61][62][63][64][65][66]. To create an event horizon, the group velocity of the modulation signals should be comparable to the propagation velocity.
In contrast to proposals which are analogs of creating cosmological particles [61,65], in this study we focus on the parameter regime where the velocity of the modulation signals is much slower than the photon's propagation speed, which was rarely considered in previous studies. The scenario is akin to realizing nonreciprocal sound propagation in an elastic metamaterial via spatiotemporally modulating the stiffness [67][68][69][70][71][72][73][74], where the Brillouin-scattering process will lead to chiral acoustic wave transport. In this work we first derive the dispersion relation of the metamaterial waveguide by employing the generalized Bloch-Floquet expansion formula. We find that there are three nontrivial dispersion regimes, which respectively correspond to left/right chiral emission and band gaps. Then we show that the metamaterial waveguide can be engineered arXiv:2206.06579v2 [quant-ph] 15 Feb 2023 as a chiral quantum channel when coupling it to superconducting qubits. The chiral transport is similar to the proposal of realizing nonreciprocal photonic transport via acoustic modulation [75]. However, the modulating signal in our proposal is the electromagnetic current bias (rather than acoustic waves), which can be tuned much faster and freely. Our proposal can be employed for demonstrating chiral quantum optics widely discussed in Refs. [76][77][78][79][80][81][82][83][84][85], and provide a novel method to realize nonreciprocal transport of microwave photons without using the bulky ferrite circulator [86,87]. Very recent studies with circuit-QED setups showed that chiral emission from qubit pairs (coupled to a common waveguide) can be realized via linear interference process [80][81][82][83].
However, the chiral transport is restricted at the single-photon level. Our proposal does not have such limitations and can maintain the chiral behaviour even in the multi-excitation regime [88][89][90][91][92]. In the last part, we show that our proposal can be extended as a chiral quantum network where remote nodes are mediated with no information back flow. The artificial waveguide is composed by a SQUID array, which impedance is modulated by a programmed current bias in each unit structure. The current signal I(j, t) changing with both SQUID's positions j and time t will spatiotemporally modulated the flux Φj through each SQUID via mutual inductance Mj. The Josephson (ground) capacitance of each SQUID is denoted as CJ (C0). The node flux φj at position j follows the Kirchhoff current relation in Eq. (2). A transmon qubit, which is in the form of a split junction with Josephson capacitance Cq, interacts with the waveguide at position x0 = 0 via a coupling capacitance C g J .
As shown in Fig 1, we consider that the Josephsonmetamaterial waveguide is composed by a SQUID array which allows microwave photons propagating along it [52][53][54]. Each node is connected to the ground with a capacitance C g . Two neighbor sites are separated with distance d 0 . The SQUID can be viewed as a nonlinear Josephson inductance L j in parallel with a capacitance C J . A series of current biases which can be programmed with external signals will produce site-dependent flux Φ j via mutual inductance M j . The relation between L j and Φ j is given by [93][94][95][96][97][98] where E s0 is the Josephson energy of one junction. Defining the node flux as φ j , we obtain the following Kirchhoff current equation for the SQUID chain [52,53] We assume that the current biases are composed of a static part and a time-dependent part, respectively. That is, where I 0 (δI) is the amplitude of the DC (AC) part satisfying δIf (j, t) where As discussed in Appendix A, given that φ j varies slowly in the length scale d 0 (i.e., long wavelength limit), the difference Equation (2) can be written in a quasicontinuous form. Defining c g = C g /d 0 (c J = C J /d 0 ) as the ground (Josephson) capacitance per unit length, the wave function in the quasi-continuous limit becomes [98] where g(x, t) is expressed as with l J being the Josephson inductance per unit length. The spatiotemporal modulation signal is encoded in f (x, t). When C J = 0, Eq. (7) is simplified as a transmission-line equation in which the impedance is modulated by a space-time-varying wave. A non-zero Josephson capacitance c J will entangle both spatial and temporal differentials together, which will produce a nonlinear dispersion for the SQUID waveguide. In  Appendix A and following discussions, we discuss its effects and give the conditions where c J can be neglected.
The eigen-wave functions of the field are obtained by adopting a generalized Bloch-Floquet expansion [68,75,99] where ω l is the quasi-energy of the lth band, and k is the wave number. For convenience, we define v 0 = 1/ l J c g as the propagation velocity of the static SQUID metamaterial waveguide without modulation. Now we consider the bias current g(x, t) varying periodically in both space and time, i.e., The simplest modulation signal is a travelling wave, and in this work we adopt f (x, t) as The phase velocity is v d = Ω d /k d with k d = 2π/λ d (Ω d = 2π/T d ) being the wave number (angular frequency). By substituting Eq. (9) into Eq. (7), the dispersion relation is obtained by solving a quadratic eigenvalue problem [68]. Detailed derivations can be found in Appendix A.
According to the experiments in Refs. [53,54,59], we list the parameters adopted for our numerical simulation in Table I. Although we can employ the generalized Floquet form in Eq. (9) to derive the spectrum, it does not mean that all kinds of modulation waves can produce stable eigenmodes in the SQUID transmission line. As discussed in Ref. [100], when the modulation velocity v d is faster than v 0 , the eigenfrequencies become complex, indicating that the oscillations are time-growing and unstable. The high-speed modulation signal can be employed for conversing vacuum fluctuations into photons, which is analog to the Hawking effect in a gravitational system [61,66].
In this study, we focus on the parameter regimes with v d v 0 . The first Brillouin zone (BZ) is limited within k ∈ (−0.5k d , 0.5k d ]. There will be Brillouin-scattering processes between modes k and k + k d (which are wave numbers for the eigenmodes of the static waveguide) due to the conservation of momentum. Additionally, the conservation of energy also requires ω l (k) ± Ω d ω l (k +k d ). Since Ω d is much smaller than ω l at k = k d /2, the modes around k ±0.5k d will interact with each other. The interactions between modes produces two band gaps, which is similar to the appearance of an anticrossing point in a two-coupled-mode system.
We assume that the waveguide is long enough (L → ∞) to support the microwave photons propagating without reflection. By quantizing the field, the root-mean-square voltage operatorV (x, t) is written as [16] V where C t = Lc g is the total capacitance of the waveguide, and a lk (a † lk ) is the annihilation (creation) operator for mode k in the lth quasi-energy band. We numerically plot ω l (k) changing with wave number k in Fig. 2(a) by adopting parameters in Table I. We find that an energy gap emerges between the two lowest bands in the first BZ (−0.5k d , 0.5k d ]. As discussed in Appendix A, under the following condition the nonlinear effect led by c J can be neglected. By adopting the parameters in Table I, one finds that k d d 0 = 0.25×10 −2 1. According to Eq. (12), even when C J = 500C g , the nonlinear effect of Josephson capacitance is not apparent, and the dispersion relation is quite close to that with C J = 0.
Equation (9) shows that the time-dependent part of φ lk (x, t) depends on both band index l and Floquet order n, and the distribution ratio of the nth Floquet order is |u ln (k)| 2 (note that n |u ln (k)| 2 = 1) for a certain mode lk. Since the modulation signal g(x, t) can be viewed as a perturbation, |u ln (k)| decreases quickly with n according to our numerical calculations. For the parameters in Table I, it is accurate enough to consider n = 0, ±1, and the corresponding dispersion relation is plotted in Fig. 2 The modulation signal g(x, t) propagates unidirectionally (rather than a standing wave carrying opposite momentums ±k d ), which will open two asymmetric energy gaps located around k = ±0.5k d (see Fig. 2). Consequently, the unconventional spectrum regime in Fig. 2 can be divided into three parts. The gray area between two quasi-energy bands, where |u ln (k)| 0 are valid for all Floquet orders {l, n}, is the conventional band gap with no propagating mode in the waveguide. In red (blue) area where point A (B) is located, the dispersion relation is asymmetric and the Floquet order {l = 1, n = 0} ({l = 2, n = −1}) possesses the highest distribution ratio |u ln (k)|. In these two regimes  9)] are shown. The distribution ratio of each Floquet order is mapped with colors. The nontrivial spectra can be divided into three parts. In the red (blue) area, the microwave photons propagate into the right (left) direction; In the grey area, the microwave photon will be localized due to no propagating modes.
microwave photons propagate unidirectionally, which will be discussed below by considering a superconducting qubit interacting with this metamaterial waveguide.

III. CHIRAL EMISSION OF SUPERCONDUCTING QUBITS
A. Interaction Hamiltonian of metamaterial waveguide coupling to superconducting qubits As depicted in Fig. 1(a), we first consider the simplest QED setup where a transmon qubit interacts with the metamaterial waveguide via a coupling capacitance C g J at x 0 . The transom qubit is composed by two identical junctions which form a SQUID loop, and its Hamiltonian is written as [12,13] wheren (φ) is the charge (phase) operator of the transmon, C Σ = C q J + C q with C q being the Josephson capacitance, and E q J (E C = e 2 /(2C Σ )) is the Josephson (charging) energy of the junction. The interaction Hamiltonian between the transmon and the waveguide is written as [13] In the limit E C E q J , the transmon can be viewed as a Duffing nonlinear oscillator. Given that only the two lowest energy levels are considered, H q can be approximately written as [13] where the charge operator iŝ In the rotating frame of ω q σ z /2 + lk ω l (k)a † lk a lk and by adopting the rotating-wave approximation, we rewrite H int as where the coupling position is set at x 0 = 0 without loss of generality, and g lk is the coupling strength which is derived as Here we take the transmon as an example. As discussed in experimental work in Refs. [59,60,98], the SQUIDmetamaterial waveguide can also interact with a flux or charge qubit, and all these circuit-QED setups can be employed to demonstrate the unconventional emission behaviors discussed in this work.

B. Chiral emission in Markovian regime
Similar to the standard spontaneous emission process, the intensity of radiation field should be narrowly centered around the atomic transition frequency ω q . Therefore, by replacing ω l (k) in Eq. (18) with ω q , g lk becomes mode-independent, i.e., g lk g 0 . We assume that a single excitation is initially in the transmon qubit, while the waveguide is in vacuum state |0 . Due to the conservation of excitation number for the Hamiltonian in Eq. (17), the system's state at time t can be expressed in the single-excitation subspace as where |1 lk represents a single photon being excited in mode a lk . The evolution governed by H int can be derived from the following coupled differential equations where ∆ l (k) = ω q − ω l (k) is the frequency detuning between the qubit and the mode a lk . By substituting the integration form of Eq. (21) into Eq. (20), we obtain where G(t, t ) is the time-delay correlation function. Given that n = n , e −i(n−n )Ω d t are fast oscillating terms, and their contributions to the evolution will be significantly suppressed when the decaying time scale is much longer than the oscillating period Ω −1 d . Under these conditions, we can only keep the resonant term n = n , which is similar to the rotating-wave approximation [75]. Finally we obtain In the emission spectrum, the intensity of the field will center at the modes satisfying the resonant condition which solutions of k are denoted as k ln . For example, in Fig. 2(b) for the frequency at the dashed position in the red regime, we mark the resonant positions (green dot A) for {l = 1, n = 0}. Around each k ln , one can approximately derive the dispersion relation as where v ln = ∂∆ l (k ln ) ∂k k ln is the group velocity. When the qubit transition frequency is far away from the band edges, G(t, t ) can be written as The delta function in Eq. (27) where Θ(x) is the Heaviside step function, Γ R(L) corresponds to the decay rate to the right (left) hand of the waveguide, and Γ 0 is the characteristic decay rate for the qubit which is derived as In the single-excitation subspace expressed in Eq. (19), the photonic wave function in real space is [2] Table I with v d = 0.05v0.
The photonic energy decaying into the right (left) hand side is expressed as By neglecting the local decoherence and dephasing of the transmon, the chiral factor is defined as [78] β As indicated in Eq. (29), the directional emission rates depend on both distributions ratios |u ln (k ln )| and group velocities v ln . For example, given that the qubit frequency lies in the red regime of Fig. 2(b), the Floquet order {l = 1, n = 0} has the largest |u ln (k ln )| (see point A), and the corresponding v ln is positive. Consequently, the photon will be chirally emitted to the right direction. In the blue regime (point B), the emission chirality will be reversed. Below we will numerically discuss the emission behaviors in different frequency regimes.

C. Numerical discussions
The chiral spontaneous emission process can well be described by the master equation, which however discards much information due to a lot of approximations. For example, master equation cannot describe the directional field distribution and Non-Markovian dynamics led by band-edge effects. To avoid those problems, we numerically simulate the unitary evolution governed by Hamiltonian in Eq. (17) by discretising the modes of the modulated waveguide in the momentum space. The photonic field in the waveguide is recovered from Eq. (31). The details about numerical methods are presented in Appendix B.
By changing the transmon's frequency, we plot the evolutions of |c e (t)| 2 in Fig. 3(a). When ω q is in the conventional dispersion regime [ω q /(2π) = 2.70 GHz], the transmon exponentially decays its energy, and the photonic field is symmetrically emitted into both left and right side [see Fig. 3(b)]. In the right chiral regime with ω q /(2π) = 3.02 GHz [see red horizontal line in Fig. 2(b)], the solutions for Eq. (25) correspond to the intersection points with the dispersion curves of different Floquet orders. The intersection point A with the Floquet order {l = 1, n = 0} has the largest distribution ratio |u ln (k)| 1, while the other |u ln (k)| are of extremely low amplitudes. The group velocity for point A is positive. Therefore, the transmon will chirally emit photons into the right part of the waveguide, as indicated by Eq. (29). Similarly, when transmon frequency is in the blue regime, most of the emitted photonic field will distribute on its left hand side. The chiral field distributions changing with time are shown in Fig. 3(c, e), respectively. Moreover, the chiral factor is about β ± 0.95, indicating that this SQUID-metamaterial waveguide can be implemented as a well-performance directional quantum bus.
When the qubit frequency lies within the band gap [the gray regime in Fig. 2(b)], the distribution ratios for all Floquet orders are around zero, i.e., |u ln (k)| 0, indicating that there is no resonant mode which can lead to exponential decay of the transmon qubit. In this scenario, the transmon only decays its partial energy into the waveguide [see Fig. 3(a) for ω q /(2π) = 3.34 GHz], and the field is localized around the coupling position without propagating outside, as shown Fig. 3(d)]. This unconventional evolution can be understood from Fig. 2(b): in the gray regime, the coefficients of all the Floquet orders are around zero, i.e., |u ln (k ln )| 0. Due to no resonant mode, only the modes which are of large detuning will contribute significantly to the dynamics. Those modes are of large density due to band edge effects, and the field will be localized in the form of bound state. The scenario is akin to an emitter being prevented from spontaneous emission when it is trapped by the band gap of a photonic crystal waveguide [101][102][103].
In this proposal the modulation signal is programmable. Therefore, the dispersion relation of the waveguide can be tailed freely, which enables us to control the qubit dynamics on demand. For example, given that the modulation phase velocity v d is switched oppositely, the chiral direction with a certain qubit frequency is also reversed. In Fig. 4(a), we plot the chiral factor β + changing with transmon frequency by setting modulation velocity as v d = 0.05v g and v d = 0.08v g , respectively. When v d = 0.05v g , the quasi-energy band l = 1 and l = 2 are separated by a finite band gap [see gray area in Fig. 2(b)]. The gap leads to the trapped bound state when the qubit frequency lies within this regime [see Fig. 3(d)]. As discussed in Ref. [68], when increasing the modulation velocity, the gap disappears, which can be found from the dispersion relation for v d = 0.08v g depicted in Fig. 4(b). The chirality will be smoothly switched from left to right without any gaps when increasing the qubit frequency. With a larger v d the chirality is enhanced and the directional bandwidth becomes wider [see Fig. 4(a)], which is due to that the Brillouin-scattering process emerges between the modes with a large energy difference. However, as discussed in Refs. [68,100], when the modulation velocity v d is comparable to v 0 , the waveguide's eigenfrequencies become complex, indicating the field is time-growing and unstable. Due to this, the modulation velocity v d should be much smaller than v 0 to avoid the unstable phenomena emerging in the whole setup.

IV. CHIRAL PHOTON FLOW BETWEEN TWO QUANTUM NODES
FIG. 5. By mediating remote nodes, the metamaterial waveguide can work as a common quantum bus in a chiral quantum network. The metamaterial waveguide is placed in a tube which is below the superconducting temperature (see experiments in in Ref. [104]). To suppress the thermal microwave noise, each node is placed in a dilution refrigerator which temperature is around T ∼ 10 mk.
By considering multiple nodes interacting with the same metamaterial waveguide, our proposal in Fig. 1 can be extended as a chiral quantum network. The schematic cascaded network is depicted in Fig. 5, where nodes are placed in separated dilution refrigerators with temperature ∼ 10 mK to suppress the thermal noise. As discussed in Ref. [104], the metamaterial waveguide can be inside a multi-sleeve tube which is below the superconducting critical temperature. This experimental method allows to connect transmons located in different cryogenic refrigerators.
Given that the transmons' frequencies are identical, the interaction Hamiltonian can be written as (34) where x i is the coupling position of the ith node. As discussed in Appendix C, we can derive the cascaded master equation for multiple nodes by tracing over the waveguide's freedoms. Taking that two transmons chirally decay/absorb the right propagating photons for example, the evolution is governed bẏ where ρ s is the reduced density matrix operator for two transmons, Γ R is the decay rate to the right side, and L = √ is the collective jump operator. The last term in H eff is unique to the cascaded quantum system, and describes the irreversible process that a photon emitted by transmon j will be reabsorbed by transmon i, while the information back flow is prevented [78]. As discussed in Appendix C, when deriving the cascaded master Eq. (35), we assume that the propagating time τ ij = (x i − x j )/v ln between two nodes is much smaller than the decaying time scale Γ −1 R . Therefore, we can adopt the Markovian approximation t − τ t, and the evolution becomes independent of time delay. This approximation is valid only when the separation distance x ij is much shorter than the wavepacket's size.
In cascaded master equation (35), the evolution information, such as the field distribution and timedecay effects, has been discarded due to taking a trace over waveguide's freedoms. To proceed beyond those limitations, our simulation is based on the unitary evolution governed by the Hamiltonian in (34), which can well describe the time-delay effects [85]. Numerical details can be found in Appendix B. In Fig. 6(a), by assuming a single excitation initially in transmon 1, we plot the evolution of two transmons with a separation distance x 2 − x 1 = 125λ d . At the frequency adopted in Fig. 3(c), the group velocity is around v ln 0.83v 0 , which is slower than v 0 due to the nonlinear dispersion relation around the band edge. Consequently the wave front arriving at transmon 2 is calculated as τ ij = 125λ d /0.83v 0 0.022ms [see Fig. 6(a)]. The field distribution |ψ γ (x, t)| 2 is depicted in Fig. 6(b), where one finds that the photon emitted by transmon 1 propagates unidirectionally, and the field absorbed by transmon 2 (inside the box) is aslo re-emitted to the right side. Due to no photonic energy back flow, transmon 1 will not be re-excited. Those numerical results show that our proposal can behave as a well-performed cascaded quantum network.

V. SUMMARY AND OUTLOOKS
In this work, we propose how to employ the Josephson array as chiral metamaterial waveguide for circuit-QED. The waveguide is in the form of SQUID chain which impedance is modulated with bias currents. When the bias signals are programmed as travelling waves, the symmetry between the left and right propagating modes is broken due to the Brillouin-scattering process. We also discuss the quantum optical phenomena by considering superconducting qubits coupling to this metamaterial waveguide. By applying the optimized modulating parameters, the qubit will emit photons unidirectionally, and chiral factor can approach 1. Last we extend our proposal as a multi-node quantum network, and demonstrate that the chiral transport between remote nodes can be realized. Compared with routing microwave photons unidirectionally with the bulky ferrite circulators, our proposal does not require strong magnetic field, and the direction can be freely tuned by the programmed bias signals.
Note that we only focus on slow travelling waves to modulate the SQUID-chain's impedance. Exploring other modulating parameters or forms, such as standingwave modulations and pulses with different shapes, might allow us to observe more intriguing QED phenomena in this metamaterial platform. As discussed in experimental studies in Refs. [56,59], the SQUID number in a single metamaterial waveguide can be around 10 3 10 4 , indicating that our proposal is within the capability of current technology. We hope that our work can inspire more efforts being devoted to exploiting SQUID metamaterials for controlling microwave photons in SQC setups. We start from the left side of Eq. (2), i.e., the difference terms related to time-dependent inductance L j (t). To derive the corresponding quasi-continuous differential form, we first rewrite it as In this work we assume that each SQUID's size is much smaller than wavelengths of both modulation wave and microwave photons. Therefore, we replace φ j (t) → φ(x, t) and f (j, t) → f (x, t). Note that G j (t) becomes The field distribution along the waveguide changes with t. At x = 125k d , the emitted photon is scattered by transmon 2. Due to directional transport, the scattered field only propagates into the right direction. The adopted parameters are the same with those in Fig. 3(c).
By replacing x = id 0 , Eq. (A1) is changed as Similarly, the left side in Eq. (2) which contains capacitance terms can also be written in quasi-continuous differential form. Finally, the nonlinear wave function of the modulated SQUID waveguide is derived as which is equivalent to the form in Eq. (7). Given that the modulation is of the travelling wave form, one can decompose the wave function in Eq. (2) as a matrix form by employing the orthogonal relations. The dispersion relation is obtained by solving the following quadratic eigenvalue problem where M 1,2 are diagonal matrices and expressed aŝ The matrixM 0 isM where We find that the Josephson capacitance c J only appears in the diagonal terms ofM 0,1,2 , from which on can easily find that c J can be neglected under the following condition Further discussions of nonlinear dispersion led by c J can be found in Sec. II of the main text.
Appendix B: Numerical methods for simulating spontaneous emission and non-Markovian dynamics Although the master equation can well describe the spontaneous emission of emitters, the Born-Markovian approximation is not valid when the interaction strength is comparable to the band width of baths. Moreover, the information of emitted photons will be traced off when deriving the master equation. To avoid those, our simulation is based on the unitary evolution governed by the original Hamiltonian in Eqs. (17,34). In the single-excitation subspace and taking for the case of two transmons for example, the system's state is written as |ψ(t) = lk c lk (t)|g, g, 1 lk +c e1 (t)|e, g, 0 +c e2 (t)|g, e, 0 , where |1 lk represents the single excitation being in lkth mode, and |0 corresponds to the waveguide in its vacuum. Since the Hamiltonian is expressed in momentum space, we first need to calculate waveguide's eigen-frequencies and wavefunctions by employing the method presented in Appendix A. Next we discretise modes in the first BZ k ∈ (−0.5k m , 0.5k m ] with a large number N , which is equal to considering a waveguide with a length L = N λ m . A similar method can be found in Ref [85]. In the single-excitation subspace, the Hamiltonian in Eq. (17) [Eq. (34)] can be mapped into a matrix with dimension N + 1 (N + 2) when the system contains a single emitter (two emitters). Taking two emitters (Q = 2) for example, the corresponding matrix is written as where φ lki (x j , t) is numerically obtained via methods in Appendix A. After obtaining the matrix form in Eq. (B1), we can numerically solve the evolution governed by the Schrödinger equation. At certain time t = t i , the amplitude c kl (t) of each mode is recorded, and the field distribution can be recovered by substituting them into Eq. (31). Employing this method, we obtained both transmon's and waveguide's evolution shown in the main text.
Appendix C: Cascaded master equation for multi-nodes system In this part we will derive the cascaded master equation for multiple emitters mediated by the metamaterial waveguide. By expanding the Schrödinger equation to the second order, the evolution of the system is expressed as [2] where ρ (ρ s ) is the density matrix operator of the whole system (the emitters), and Tr R represents taking a trace over the waveguide's freedoms. We assume that the waveguide is always approximately in its vacuum state. Therefore, the correlation relations for the waveguide's modes satisfy a lk = a † lk = a † lk a lk = 0 and a lk a † lk = 1. By substituting those relations into Eq. (C2), we obtainρ where the correlation function A (x i , x j , t, t ) is defined as [75] A (x i , x j , t, t ) = e iωq(t−t ) lk φ * lk (x i , t)φ lk (x j , t) with x ij = x i − x j being the distance between two emitters. We have neglected the phase terms e i(k ln +k d )(xi−xj ) in A (x i , x j , t, t ). Note that the integral regime is within t ∈ [t 0 , t]. Therefore, the δ-funtion will produce non-zero value in Eq. (C3) only under the condition t − x ij /v ln ≤ t. The delay-time τ ij = x ij /v ln (i = j) corresponds to the propagating time between two separated emitters. When x i > x j (x i < x j ), the nonzero time-delay correlation is mediated by the right (left) propagating modes with v ln > 0 (v ln < 0). Given that τ ij is much smaller than the time scale of spontaneous emission, we can replace t − τ → t, and Eq. (C3) is simplified aṡ where Γ R(L) is the decay rate into right/left side, which are expressed in Eq. (29). In Eq. (C5) we have employed the properties of δ-function When the chiral factor of the whole system approaches β ± 1, we obtain the cascaded master equation in the main text.