Supplemental Materials: Frustration-induced anomalous transport and strong photon decay in waveguide QED

Ron Belyansky, 2 Seth Whitsitt, 2 Rex Lundgren, 2 Yidan Wang, Andrei Vrajitoarea, Andrew A. Houck, and Alexey V. Gorshkov 2 Joint Center for Quantum Information and Computer Science, NIST/University of Maryland, College Park, Maryland 20742 USA Joint Quantum Institute, NIST/University of Maryland, College Park, Maryland 20742 USA Department of Electrical Engineering, Princeton University, Princeton, NJ, USA (Dated: February 19, 2021)

I. Renormalized spin frequency 1 II. Numerical methods 2 S1. Orthogonal polynomial mapping 3 III. Elastic S-matrix in terms of spin susceptibilities 4 IV. Derivation of the spin susceptibilities 5 V. Inelastic scattering 9 VI. Comparison with a model without frustration 11 References 13

I. RENORMALIZED SPIN FREQUENCY
In this section, we reproduce the derivation [S1] of the crossover energy scale between weak and strong coupling, the equivalent of the Kondo temperature, which we associate with the renormalized spin frequency ∆ R in Eq. (2) in the main text.
We employ the renormalization group (RG) procedure as used in Ref. [S1]. Let us denote the cutoff at some energy scale l by Λ(l) (such that Λ(0) ≡ ω c is the original cutoff of the problem). The RG procedure consists of integrating out the high-energy modes, and thus redefining the cutoff from Λ to Λ − dΛ. This leads to the RG flow equations for the coupling constants α and h ≡ ∆/Λ [S1] (that can also be derived from perturbation theory in α) where dl = −dΛ/Λ, which implies that l = log ωc Λ . The equation for α can be readily solved, (S3) α(l) = α(0) 1 + 2lα(0) .
Note that both the cutoff Λ and ∆ are decreasing as function of l. Equivalently, h increases from its initial value of h(0) = ∆/ω c . Eventually, we have h(l * ) ≈ 1, which occurs when ∆(l * ) ≡ ∆ R = Λ(l * ), and hence the RG breaks down. The energy scale corresponding to this is l * = log ω c /∆ R . Plugging this into Eq. (S4) gives Eq. (2) from the main text, . (S5) This differs from the result presented in Refs. [S1, S2] in that we kept the ∆ R on the right-hand-side (whereas these references approximated it by ∆) as it agrees better with the numerical results.
As we show later in Sec. IV, we also reproduce exactly this equation by applying the Callan-Symanzik equation to the bare Green's function of the spin.

II. NUMERICAL METHODS
In this section, we describe the numerical methods we use to compute the single photon transport properties. We first write the Hamiltonian (Eq.1 of the main text) in terms of bosonic creation and annihilation operators, as follows, . We make the transformation which transforms the Hamiltonian into two commuting partŝ whereg(k) = √ 2g(k). It is enough therefore to only simulate the dynamics of theĤ XY SB Hamiltonian. Explicitly, to determine the single-particle scattering properties, we perform the following procedure. We create a single-particle Gaussian wavepacket on top of the ground state, with amplitude c k = N e − (k−k 0 ) 2 2σ 2 +ikx0 , where N is a normalization so that kmax 0 dk|c k | 2 = 1. Without loss of generality, we choose this excitation to be of the x type (since for α x = α y the Hamiltonian is invariant under x ↔ y exchange.). We then evolve this state in time, which leads to the scattering of the wavepacket off the spin. At long times after the scattering, we can extract several observables such as the elastic scattering amplitudes and number of elastic and inelastic photons in the final state.
Denoting the ground state of the full HamiltonianĤ by |GS = GS XY SB ⊗ 0 f ree , the initial state is The time-evolved state is Thus we see that all quantities of interest can be obtained from correlation functions and matrix elements of the states GS XY SB and |ψ XY SB (t ∞ ) .

S1. Orthogonal polynomial mapping
The HamiltonianĤ XY SB from Eq. (S8) describes a system with very nonlocal interactions. In order to efficiently simulate it with matrix-product-states, we use the orthogonal polynomial mapping introduced in [S4, S5], which maps the Hamiltonian into a tight-binding model with only nearest-neighbor interactions.
Here we summarize the main steps of the mapping. For more details, see Refs. [S4, S5]. We work with the Hamiltonian from Eq. (S8), reproduced here (S16) where ω(k) = ω c k, k max = 1, andg i (k) = √ 2α i ω c k. The resulting spectral functions are We introduce the unitary transformation where p i,n (k) are orthonormal polynomials with respect to the measure dµ i (k) =g 2 i (k)dk [i.e. p i,n , p j,m ≡ kmax 0 dµ i (k)p i,n (k)p j,m (k) = δ nm δ ij ], and a set of new discrete bosonic modes that satisfy b i,n ,b † j,m = δ ij δ n,m . Using the recurrence relations of orthogonal polynomials, one can show that the Hamiltonian in Eq. (S16) can be written as which describes two semi-infinite tight-binding bosonic chains that are both coupled to the spin via their first site. For the Ohmic spectral function, the p i,n polynomials are the Jacobi polynomials, and the on-site energies and hopping coefficients are given by ν n = 1 2 1 + 1 (1 + 2n)(3 + 2n) , β n+1 = 1 + n 1 + n + 2n 1 + n 2 + n . (S21) we can convert measurements in theb basis to observables in frequency space. For example, the frequency-mode occupation is (for k = k ) Note that this is an exact mapping, provided the length of the chains is infinite. In practice, the length of the chains is truncated to a finite value L, and the dimension of each bosonic Hilbert space is truncated to a finite value d.
We varied these parameters and found that L = 250, d = 5, and bond dimension of χ = 30 are adequate to obtain converging results for the scattering for most values of α.

III. ELASTIC S-MATRIX IN TERMS OF SPIN SUSCEPTIBILITIES
In this section, we derive the relation between the elastic scattering coefficients and the spin susceptibilities, given in Eqs. (3,4) in the main text.
Let us write the Hamiltonian aŝ We are interested in the S-matrix between a particle with momentum k in bath i and a particle with momentum k in bath j: where ψ ± ik are the exact incoming and outgoing scattering eigenstates. Following Ref. [S6], we can write these eigenstates as follows: where |ψ 0 is the ground state of the fullĤ with energy E 0 , and χ ± i are the states of the scattered particles. Since ψ ± ik are eigenstates ofĤ with energy E 0 + ω k , Schrodinger's equation implies Substituting this into Eq. (S29) gives with > 0 taken to zero at the very end. From this we find Plugging this back into the S-matrix Eq. (S27) gives where we defined the T -matrix To evaluate the first term, we perform a similar manipulation as in Eq. (S30): Inserting this into Eq. (S34) yields which, when ω k = ω k , we recognize as the Fourier transform of the retarded Green's function which we equivalently refer to as the spin susceptibility in the main text.

IV. DERIVATION OF THE SPIN SUSCEPTIBILITIES
In this section, we explicitly compute the spin susceptibility, Eq. (S39). We will do so by first computing the imaginary-time Matsubara Green's function, Here, the imaginary-time dependence of operators isσ j (τ ) = e Hτσ j e −Hτ , and T τ is the imaginary-time ordering operator. This function may be computed using the usual diagrammatic perturbation theory if we can use Wick's theorem, after which we may obtain the spin susceptibility by analytic continuation [S7]: However, the Pauli matrices do not satisfy Wick's theorem. We can get around this by using an Abrikosov pseudofermion representation of the spins. We introduce a two-component set of fermions, {χ a , χ † b } = δ ab , a, b = 1, 2, related to the spin operators byσ This is only a faithful representation of the spin operators in the subspace χ † 1 χ 1 + χ † 2 χ 2 = 1. However, we can project to this subspace using the Popov-Fedotov trick of using an imaginary chemical potential µ = −iπ/2β, which results in a cancellation between the unphysical subspaces [S7]. This method requires us to work at finite temperature during intermediate calculations, but below we will always take the β → ∞ limit as early as possible. We may now express the system as a coherent-state path integral. We introduce the bosonic fields φ ±,k = (a x,k ± ia y,k )/ √ 2, after which our model may be described by the Lagrangian In this form, it is straightforward to treat the interactions g k perturbatively using a Feynman-diagram expansion and the Matsubara formalism. We have the bare fermionic Green's functions with where α n = π(2n + 1)/β, n ∈ Z. Similarly, the bosonic propagators are and Ω n = 2πn/β, n ∈ Z. The interaction terms in Eq. (S43) result in four interaction vertices, which contribute dependence on g k . Since we are interested in correlation functions of the spin operators, we also introduce a diagrammatic notation representing the composite operatorsσ + = χ † 1 χ 2 andσ − = χ † 2 χ 1 . The Feynman rules for this theory are displayed in Fig. S1.
As a demonstration of this formalism, we obtain the spin susceptibility in the non-interacting (α = 0) limit by computing the diagrams with a single fermion loop and no bosonic propagators, We also have G −+ (iΩ n ) = G +− (−iΩ n ) and G ++ = G −− = 0. Going back to the xy basis, taking β = ∞, and analytically continuing, we obtain the expected form for the spin susceptibilities for α = 0: These expressions could be simply obtained through a direct computation at zero temperature with the spin operators, but the diagrammatic approach is useful for including interactions. We note that the susceptibilities in Eq. (S49) have a simple pole located at the bare spin frequency, and such a pole can never be shifted or broadened by computing a finite number of diagrams. Therefore, we will use both the Callan-Symanzik equations and a Dyson equation to sum an infinite number of diagrams, which will result in a change in the analytic structures of the susceptibilities.
As discussed in Sec. I, if we perform an RG transformation on our system, redefining the cutoff from ω c ≡ Λ(0) to some lower cutoff Λ(l), we obtain the flow equations where h ≡ ∆/Λ. In addition to coupling renormalization, it turns out that one also needs to renormalize the spin operators under an RG transformation, and one may show that, perturbatively, implying a flow for the operators, These flow equations, first obtained in Ref.
[S1] using a Wilsonian momentum-shell RG scheme, may also be obtained by treating Eq. (S43) using the conventional methods of quantum field theory. We now use the fact that the susceptibilities should be independent of an RG transformation, dG R ji /dΛ = 0. Taking into account any explicit and implicit dependence on the cutoff, this implies the Callan-Symanzik equation, The general solution to this partial differential equation is where the f ji are arbitrary functions of the "running couplings," defined as Comparing this general solution to the leading-order expressions of Eq. (S49), we may read off the α = 0 limit of the functions f ji , and then use the α-dependence implied by the solution of the differential equation to find where we have plugged in Λ(l = 0) = ω c to give expressions in terms of the initial cutoff and the bare quantities, and ω has a small positive imaginary part. We see that both expressions no longer diverge at ω = ∆, but instead they have poles at ω = ∆ R , where ∆ R satisfies From Eq. (S57), we see that the effect of solving the Callan-Symanzik equation was to sum the "leading logarithms," which are terms of the form α n log n ω c /ω at nth order in perturbation theory. Although we have succeeded in capturing the renormalization of the spin frequency using the Callan-Symanzik equations, they still predict a sharp behavior at ω = ∆ R , whereas we expect interactions to broaden the peak near the renormalized spin frequency. We rectify this by computing an additional infinite set of diagrams using Dyson's equation. We first consider the one-particle irreducible Green's functions, G 1PI ji , which are defined to be the complete set of diagrams that cannot be split in two by cutting a single propagator. By the structure of the interactions, the only possible propagators which can be cut to disconnect a susceptibility diagram is a bosonic propagator. As a result, we have the exact relation (in matrix notation) Then the full Matsubara Green's function satisfies We now approximate the full Green's function by just using the leading-order result, Eq. (S48), for G 1PI . This corresponds to summing all possible "bubble diagrams" which contribute to the susceptibility, which is reminiscent of the RPA approximation in the dense electron gas. In this approximation, we obtain the susceptibilities as We see that the inclusion of these diagrams has resulted in a finite imaginary part in the denominator, which removes the pole on the real-ω axis. We may now furthermore use the Callan-Symanzik equation to sum the leading logarithms. After matching Eq. (S55) to Eq. (S62), we obtain the spin susceptibilities given in Eqs. (5,6) of the main text. As noted in the main text, to this order, we have found that the above expressions do not lead to any contribution to inelastic scattering (γ(ω) in the main text). We have checked that including all O(α) contributions to G 1PI still does not lead to inelastic scattering. We believe that including O(α 2 ) contributions to G 1PI will lead to inelastic contributions, which is consistent with Sec. V, where we show that inelastic contributions to the S-matrix appear at O(α 2 ).

V. INELASTIC SCATTERING
In this section, we will consider the leading contributions to inelastic scattering in perturbation theory using the diagrammatic approach developed in Sec. IV. To proceed, we need a relation between time-ordered expectation values and S matrix elements. Such a relation is called the LSZ reduction formula in relativistic quantum field theory [S8], but we can follow the derivation for our present system and derive a non-relativistic analogue of the reduction formula. If we consider the scattering of n photons with momenta k 1 , k 2 , ..., k n into a state with n photons with momenta k 1 , k 2 , ..., k n , the S matrix element is given by Here, ω k = |k| is the energy of the photon. Note that the real time-ordered correlation function appears in this expression, which is related to the Matsubara correlation functions in Sec. IV by a Wick rotation. This expression greatly simplifies after Fourier transforming to frequency space. When we evaluate diagrams using Wick's theorem, we will come across the following bosonic contractions from the external legs of Feynman diagrams, These will set the external legs of the Feynman diagrams to their on-shell values, ω = ω k , and additionally cancel out the contribution of the external propagators. Then the S matrix elements simply become That is, we evaluate the diagram in momentum space with on-shell external legs and omit the external propagators (i.e. we "amputate" the legs).
The symmetry of the model greatly reduces the number of diagrams we need to consider. In particular, the Hamiltonian has a U(1) symmetry [S9], and the conserved charge [which is more easily apparent in terms of the ± photonsâ ± = 1 2 (â x ± iâ y )] is FIG. S3. The nonzero diagrams contributing to the scattering of one photon into three photons.
Since the scattering cannot change the spin (as this would require the existence of bound spin-photon eigenstates), it follows that the quantity k (â † −,kâ −,k −â † +,kâ +,k ) must be conserved in any scattering process. We now consider the simplest case, where one photon scatters into three photons. The discussion above means that there are only two nonzero diagrams, pictured in Fig. S3. An explicit calculation finds that to lowest order we have together with the energy conservation condition ω = ω 1 + ω 2 + ω 3 (enforced by a delta function, which has been omitted in these expressions). Denoting Eq. (S67) by a function f (ω 1 , ω 2 , ω 3 ; ω), we see that Eq. (S68) is simply f (−ω 1 , −ω 2 , −ω 3 ; −ω). Alternatively, each amplitude can be found from the other by substituting ∆ → −∆ (up to an overall minus sign). This is because the two diagrams in Fig. S3 are related to each other by replacing + photons with − photons and − photons by + photons. From Sec. IV we see that this requires flipping the sign of ∆. Converting Eqs. (S67) and (S68) to the xy basis, we find the four amplitudes The resulting probabilities are: for x → xxx γ xxx (ω 1 , ω 2 , ω 3 ; ω) = α 4 ω 4 This expression is equivalent to the leading-order result given in Ref. [S10] for the spin-boson model. This can be understood from the fact that, at leading order, the x → xxx process does not involve any y photons. Computing the other scattering probabilities, we find that the processes x → {yyy, xxy} have a simple relation to the above process, given by The remaining process, x → xyy, does not have a simple relation to the above expressions. It is explicitly given by As with the elastic probabilities, the leading-order calculation leads to poles in the scattering probabilities at the bare spin frequency ∆. We may apply procedures like those in Sec. IV to obtain an analytic dependence, which better resembles that of the fully-interacting problem. For example, by applying the Callan-Symanzik equation to the amplitudes in Eq. (S67)-(S68), we will find that the instances of ∆ will be corrected to the renormalized value ∆ R . Similarly, if we developed a Dyson equation for this amplitude, we expect that the poles will be softened to broad peaks with a similar width to the peaks seen in the elastic probabilities given in the main text.
From Eq. (S74), we see that scattering processes involving a photon from one waveguide into three photons in the other waveguide will dominate over scattering entirely within the same waveguide if ω ∆ R . In this same limit, we do not expect a large region of phase space with very large final ω 3 , so γ xxy is expected to be much smaller than γ yyy . We have verified by numerically integrating the above expressions over the possible final frequencies that, when ω ∆ R , the total cross section for the processes x → {xyy, xxy} are of the same order of magnitude, and they are both much smaller than x → yyy. We also found that the cross section for x → xxx is much smaller than the three other processes in the same limit. This is consistent with our numerical results, where we found that the inelastic scattering for ω ∆ R is dominated by scattering from one photon flavor to the other.

VI. COMPARISON WITH A MODEL WITHOUT FRUSTRATION
In this section, we compare the results of the main text to the more common situation in waveguide QED, without frustration. We assume the same geometry as in the main text [shown in Fig. S4], with the spin coupling to two waveguides, but the coupling is via the same operatorσ x . The coupling constant to each waveguide is given by S4. Schematic of the same model as in the main text but without frustration. Here the spin-1/2 is coupled locally to two independent waveguides with the same operatorσx.
We obtain the elastic scattering probabilities shown in Fig. S5 by the same method as described in the main text and in Secs. II and IV above. For the analytical calculation, the only major difference is that the φ + and φ − bosons are now equivalent, so the bottom two vertices in Fig. S1(b) are equivalent to the top two. The elastic scattering coefficients shown in Fig. S4 take the form t 11 (ω) = 1 + r 11 (ω), and the susceptibility computed using the methods of Sec. IV is given by (S79) Figure S5 shows that the elastic response has a resonance that is in excellent agreement with the well-known result from the spin-boson [S11] or Kondo literature: However, we see that, in contrast to the frustrated case, the transmission away from the resonance is very high, whereas the reflection is only nonzero around the ∆ R resonance. Note that Im(χ xx (ω)) in Eq. (S79), describing the spectral weight of the spin, shows a sharp peak at ∆ R . Moreover, the large ω limit of Eq. (S79) is Re(χ xx (ω)) ∼ ω α T −2 and Im(χ xx (ω)) ∼ ω 2α T −3 , both of which decay to zero much faster than in the frustrated case. These results are consistent with previous studies [S10, S12-S15], in other one-dimensional realizations of the spin-boson model. In all of these cases, the system is well described in terms of the polaron or dressed-spin picture. Even at very large couplings, where ∆ R → 0 and the resonance is disappearing, the photons show almost no trace of the coupling to the spin, being almost fully transmitted. This should be contrasted to the frustrated case discussed in the main text, where we have the complete opposite scenario, where the spin becomes extremely widespread over the entire energy spectrum, leading to strong elastic response in the whole range ω > ∆ R .
Next, we look at the inelastic scattering, employing the same numerical procedure as we used in the main text. We scatter narrow wavepackets and record the resulting number of elastic and inelastic particles, shown in Fig. S6. Note that because the coupling to both waveguides is the same, the number of inelastic photon emitted in each of the two waveguides is also the same. Hence, we present only the total number of elastic and inelastic photons. We see that the number of elastic photons is always near 1, never going below ∼ 0.9 for the wavepackets considered. In fact, since ∆ R /∆ ≈ 0.5 for α T = 0.25 (see Fig. S5), most of Fig. S6 is in the regime ω ∆ R . The number of inelastic particles does not exceed ∼ 0.36, which occurs for the lowest-energy wavepacketω in = 0.5∆. In fact, as shown also in Ref. [S10], for a given α T , the inelastic scattering rate peaks at an energy close to ∆ R . For example, for α T = 0.5 (the Toulouse point), the peak occurs at 2∆ R [S10]. This should be contrasted again to the situation with frustration, described in the main text, where we found that the inelastic rate remained saturated close to its maximum allowed value 0.5 even for energies above the bar spin gap ∆.