Inelastic decay from integrability

A hallmark of integrable systems is the purely elastic scattering of their excitations. Such systems possess an extensive number of locally conserved charges, leading to the conservation of the number of scattered excitations, as well as their set of individual momenta. In this work, we show that inelastic decay can nevertheless be observed in circuit QED realizations of integrable boundary models. We consider the scattering of microwave photons off impurities in superconducting circuits implementing the boundary sine-Gordon and Kondo models, which are both integrable. We show that not only inelastic decay is possible for the microwave photons, in spite of integrability, and thanks to a nonlinear relation between them and the elastically-scattered excitations, but also that integrability in fact provides powerful analytical tools allowing to obtain exact expressions for response functions describing the inelastic decay. Using the framework of form factors, we calculate the total inelastic decay rate and elastic phase shift of the microwave photons, extracted from a 2-point response function. We then go beyond linear response and obtain the exact energy-resolved inelastic decay spectrum, using a novel method to evaluate form factor expansions of 3-point response functions, which could prove useful in other applications of integrable quantum field theories. We relate our results to several recent photon splitting experiments, and in particular to recent experimental data that provides evidence for the elusive Schmid-Bulgadaev dissipative quantum phase transition.


I. INTRODUCTION
Integrability entails exceptional consequences.An extensive number of local conservation laws, the defining feature of an integrable system, are at the heart of the celebrated Bethe ansatz [1], which allows for rare exact solutions of 1-dimensional interacting many-body quantum systems.Its fundamental ingredient is a set of elementary excitations whose scattering is purely elastic -in any scattering process, the number of excitations is conserved, as well as the set of their individual energies and momenta.This striking feature has drawn large theoretical interest ever since Bethe's seminal work, promoting extensions of the method which apply to a large variety of discrete and continuous 1-dimensional models [2][3][4][5][6][7][8][9][10][11].
Remarkably, integrability is no longer just a theoretical curiosity.Recent advances in the fabrication techniques of quantum simulators have enabled the experimental realization of integrable systems, leading to an interplay between experiment and theory.The past decade has seen several theoretical breakthroughs concerning the equilibrium and out-of-equilibrium dynamics of integrable systems [12][13][14][15] that go hand-in-hand with surprising experimental observations [16][17][18][19].Essentially, integrability gives rise to counterintuitive experimental measurements which push the boundaries of well-established theoretical frameworks, and improve our understanding of the role of integrability in an ever-growing list of mechanisms.
Experiments on the quantum simulation of many-body quantum models, both integrable and non-integrable, have been mostly restricted to the realm of cold atom systems.Another possible platform for quantum simulation * burshtein2@mail.tau.ac.il is that of superconducting circuits.The rapidly evolving field of circuit quantum electrodynamics (cQED) deals with the simulation of interacting models by means of Josephson junctions or their flux-tunable counterparts, the SQUIDs.Yet experimental realizations have up to now been quite limited, and mostly dedicated the use of superconducting circuits to the implementation of nonlinear bulk models exhibiting non-ergodic behavior [20][21][22][23][24][25][26][27].The field of cQED reveals its true strength in the simulation of quantum impurity models.The intrinsically large kinetic inductance of Josephson junctions allows one to design transmission lines with impedances on the order of the resistance quantum [28][29][30], providing an environment for photons with an effective fine structure constant of order unity.Furthermore, the nonlinearity of the junctions provides the means to realize many types of quantum impurities [31][32][33][34] which are strongly-coupled to the photonic environment.Single-photon spectroscopy then provides highly sensitive tools to investigate the fine details of the boundary models of interest, across a wide range of parameters, and probe fundamental phenomena in those many-body systems.
The starting point of this work is an experimentally observed phenomenon that is seemingly at odds with integrability.Recent experiments have demonstrated that photons propagating in a high-impedance transmission line, setting an environment with a large effective lightmatter coupling, scatter inelastically off a quantum impurity with a very high probability [33,34].These observations were reproduced in Ref. [35].It appears that such photon splitting, depicted in Fig. 1, has nothing to do with integrability; as stated above, scattering processes in integrable systems are highly restricted by the extensive number of local conservation laws, forbidding particle production.The flexibility of the circuit elements provides us with tools to check this assumption, HI ω γ(ω ′ |ω) γ(ω) Figure 1.(a) An incoming photon with frequency ω, injected from the antenna on the left and propagating through the transmission line in the center, may decay inelastically as it scatters off the impurity on the right, in spite of the purely elastic reflection of the fundamental excitations of the integrable models.Each photon may be represented as a combination of eigenstates composed of excitations | ⃗ λ⟩ ϵ λ with types ϵ λ and rapidities ⃗ λ, where the weights are determined by the form factors fϵ λ ( ⃗ λ), and the total energy of the excitations in each eigenstate, ν λ ∼ i e λ i , is equal to the photon frequency ω (see Section III for notations and definitions).The excitations scatter elastically off the boundary (quantum impurity), picking up phases determined by the reflection matrix R ϵ ′ λ ϵ λ ( ⃗ λ), such that the outgoing combination no longer represents a single-photon state, but rather a multi-photon state.The measured observables -the total inelastic decay rate γ (ω) and the energy-resolved inelastic decay spectrum γ (ω ′ |ω), as well as the elastic phase shift δ (ω) (not depicted here) -all shed light on the fundamental properties of the impurity models.(b) Implementation of Eq. ( 1) with the bsG and Kondo impurities (Eqs.( 3) and (5), respectively) in a cQED setup.The array of Josephson junctions and capacitors implements a high-impedance transmission line, thanks to the kinetic inductance of the Josephson junctions.
by tuning the impurity parameters to those of integrable boundary models.Strikingly, experiments show that inelastic decay persists even when the impurity parameters are pushed towards those of the boundary sine-Gordon (bsG) model [36], which is known to be integrable [37].This appears to contradict the defining feature of integrability, and raises a fundamental question -how can inelastic decay emerge in a system governed by purely elastic scattering rules?
Previous theoretical treatments of such scattering experiments used different forms of weak/strong coupling expansions [31,32,35,[38][39][40][41][42] which, despite providing successful quantitative predictions, are intrinsically limited in two manners.First and foremost, weak or strong coupling expansions do not address the apparent discrepancy between integrability and inelastic decay; in order to settle it, one needs to explicitly use the framework of integrable systems, and show how the purely elastic excitations underlying the theory can give rise to the observed photon splitting.We note that inelastic effects have also been considered both theoretically and experimentally in integrable strongly correlated electronic quantum impurity systems [43][44][45], yet a direct consideration of the elastically scattered excitations from the integrability picture has been missing.Second, in case the impurity term in the Hamiltonian is relevant, perturbation theory and strong coupling expansions are only applicable either above or below a certain renormalization group (RG) energy scale, regardless of the impurity's strength.Wilson's numerical RG (NRG) can in principle cover the entire frequency range, but is not easy to apply accurately for bosonic baths [46].This again calls for the use of the integrability framework, as it provides powerful analytical tools allowing for an exact solution at all frequencies, linking the scaling laws above and be-low the RG scale.An exact low energy solution is especially desirable for the bsG model in the context of the Schmid-Bulgadaev transition [47,48] -a 40-yearsold predicted quantum phase transition, whose lack of clear experimental proof sparked a recent debate [49][50][51][52][53][54][55][56].A low energy theory of the scattering of photons in the bsG model could supplement a recent experimental study, which seeks signatures of the transition in the inelastic and elastic scattering rates of microwave photons [36].
In this work, we show how inelastic decay can be described via the language of integrability.The principle idea is summarized in Fig. 1(a), which depicts a generic integrable quantum field theory realized in a cQED setup.We rely on the nonlinear relation between the microwave photons, which are observed to scatter inelastically, and the elementary excitations of the integrable system, whose scattering is purely elastic.This nonlinear relation is encoded in the form factors of the models -the matrix elements of the bosonic field operator in the basis of the fundamental excitations of the integrable theory [57].Building upon Refs.[58][59][60], we use the form factors to obtain exact expressions for the reflection coefficient of the microwave photons, which encodes the inelastic and elastic decay rates.We then go beyond linear response and calculate the exact energy-resolved inelastic decay spectrum, by devising a method to evaluate a 3-point response function using form factors.This is a non-perturbative and rapidly-convergent diagrammatic approach with clear physical intuition, analogous to Wick's theorem in free theories.Our technique yields more general results as compared to existing methods for the calculation of form factor expansions of multipoint correlation functions [61][62][63], and could prove useful in other contexts of integrable quantum field theories.
We apply our methods for the bsG and Kondo models, both of which qualify as special cases of Fig. 1(a), with implementations depicted in Fig. 1(b).We illuminate the fundamental physics governing each model through the lens of the exact scattering rates, and show how they act as probes which discriminate between the two models, as well as between the regimes of low and high energies.
In particular, we analyze our non-perturbative results in light of the Schmid-Bulgadaev transition, and show how one may observe signatures of the transition in the calculated rates, which have been recently measured in an experiment [36].
The rest of the paper is organized as follows.We begin by describing the bsG and Kondo models in Section II, discussing their realizations in a cQED environment and the observables commonly measured in experiments.In Section III, we summarize the quintessential features and tools used in the treatment of massless integrable quantum field theories with a boundary, and introduce some notations that are utilized in this work.Armed with this analytical power, we calculate the exact reflection coefficient of the microwave photons in Section IV, from which we extract the total inelastic decay rate and elastic phase shift.We show how inelastic decay can emerge from the purely elastic scattering rules of the fundamental excitations, discuss the key features and asymptotic behavior of our results, and relate them to recent experiments.The exact energy-resolved inelastic decay spectrum is obtained in Section V, where we introduce a method to calculate a 3-point response function using form factors.The technical details of our method are discussed in Section V B and the Appendixes, the physical intuition underlying our diagrammatic approach is emphasized in Section V C, and the advantage of the decay spectrum over the total inelastic decay rate as a diagnostic tool is explained in Section V D. We conclude in Section VI.

II. MODELS AND OBSERVABLES IN A CQED SETUP
In this work, we consider two different yet related quantum impurity models -the bsG and Kondo models.We begin by showing how they can be derived from cQED setups.The implementation of these models using superconducting circuits allows one to probe their fundamental properties in experiments via spectroscopy, which measures the elastic and inelastic scattering rates of microwave photons.
Consider a very long array of N +1 ≫ 1 superconducting grains with lattice spacing a (total length ℓ = N a), linked to each other by Josephson junctions E line J and capacitors C line , and to the ground by capacitors C g .The circuit Hamiltonian reads (setting ℏ = e = 1) where φn and Qn are the superconducting phases and charges of the grains, respectively, HI is the Hamiltonian of the impurity at n = 0, and C −1 is the inverse of the capacitance matrix C, which is given by ) in the bulk.The generic circuit implementing Eq. ( 1) is depicted in Fig. 1.Choosing the line parameters such that , where E line C = 1/ 2C line , both phase slips and anharmonic effects are strongly suppressed in the intergrain Josephson junctions, leading to an approximately quadratic transmission line.In the thermodynamic and continuum limits, a → 0 and N → ∞ with ℓ = N a fixed, the Hamiltonian reads (2) where ρ (x) = Qn=x/a /a is the charge density field.The array velocity is given by v = a/ L line C g , and is the effective inductance of the array junctions.The Luttinger parameter of the system is given by the dimensionless line impedance, is the resistance quantum.The kinetic inductance of the Josephson arrays, L line , enables the implementation of large impedances, z ∼ 1, two orders of magnitude larger than typical impedances in classical LC resonator lines [28,64].The impedance plays the role of an effective fine structure constant, determining the strength of lightmatter interaction in such cQED setups; thus, achieving z ∼ 1 is the key to high-probability inelastic decay on the single-photon level [33,34].

A. The boundary sine-Gordon Hamiltonian
The bsG model can be implemented by connecting the transmission line to a Cooper pair box, where the charging energy , which is the elastic linewidth of the impurity, equal to one over the RC time defined by the line impedance and impurity capacitance.Another high-energy cutoff is imposed by the plasma frequency of the Josephson array, ω p = 1/ √ L line C line ; the UV cutoff is then set by the smaller of these two scales, Λ ∼ min {Γ 0 , ω p }, and is assumed to be much larger than any other energy scale in the problem.We consider the scaling limit, E J , Λ → ∞, such that the ratio J is an emergent RG scale, marking the characteristic frequency below which the impurity term cannot treated as a perturbation, and the phase at the boundary φ (x = 0) is pinned to one of the minima of the boundary cosine (note that this implies that perturbation theory remains valid at all frequencies for z > 1).At zero temperature, E ⋆ J is the only remaining energy scale in the problem.In this limit, the capacitive term is effectively eliminated from the impurity Hamiltonian.Redefining the fields as ϕ = φ/ √ 2z, ρ = ρ√ 2z, we find Taking the limit ℓ → ∞, we arrive at the Hamiltonian of the bsG model.

B. The Kondo Hamiltonian
There are several ways to arrive at the Kondo model from cQED Hamiltonians.One is to connect the transmission line to a fluxonium [65] at external half flux quantum, with the impurity Hamiltonian This Hamiltonian realizes a double-well potential for the flux φf .In Appendix A, it is shown that the line and the fluxonium give rise to the spin-boson Hamiltonian, where S z is a pseudo-spin operator whose two eigenvalues, ±1/2, correspond to the two potential wells.The S x term describes the tunneling between the wells, where J is the tunneling matrix element, equal to the splitting between the two lowest energy levels of the fluxonium.The S z term describes an effective capacitive coupling of the fluxonium to the array.Note that the coupling coefficient is written in terms of some parameter z ′ , which is proportional but not identical to the normalized impedance z (see Appendix A for details); for brevity, from here on we revert to writing z instead of z ′ .The unitary transformation H → U −1 HU with U = e iϕ(x=0,t)Sz leads to the Kondo Hamiltonian [66], A similar implementation, using a flux qubit instead of a fluxonium, has been realized in Ref. [67].The same Hamiltonian can also be derived from the effective lowenergy description of two transmission lines coupled by a Cooper pair box, with Josephson junctions connecting the capacitors to the lines [31], although this implementation could suffer from some drawbacks [68].A slightly different transformation, V = e √ 2ziϕ(x=0,t)Sz , eliminates the capacitive coupling term: x=0) .(8) This will be the form used in the following, again in the limit of a semi-infinite lead, ℓ → ∞.Introducing a cutoff frequency Λ to our model, the scaling limit reads again

C. Definition of the scattering rates
Both Hamiltonians may be written in the form The bulk Hamiltonian is quadratic and diagonalized by plane waves, with the dispersion relation ω q = vq and mode spacing ∆ = πv/ℓ, so that q = n∆ where n > 0 is an integer.The mode expansion of the fields reads where b † q , b q are bosonic creation and annihilation operators.
The implementation of the Hamiltonians ( 4) and (8) in cQED setups provides direct access to their properties in a controlled environment.In this work, we focus on scattering experiments at zero temperature, where a single microwave photon at frequency ω is injected from the open end of the system (the antenna on the left in Fig. 1), and propagates towards the impurity.One may then measure response functions; the simplest one is the reflection coefficient r (x, x ′ ; ω) of a photon with frequency ω, injected at x ′ and measured at x.The dependence on x, x ′ should amount to a trivial phase term, r (x, x ′ ; ω) = e −iω(x+x ′ )/v r (ω), and r (ω) defines the inelastic decay rate γ (ω) and elastic phase shift δ (ω): There are several ways to calculate the reflection coefficient.One is to solve a scattering problem for the microwave photons, and read off r (ω) from the scattering matrix of the theory [69].Here we adopt a different approach, and calculate r (ω) from the conductance in the half-infinite line, by means of a Landauer-Buttiker-like formula.The AC conductance is given by the Kubo formula: Introducing the right and left current operators, we find The terms G RR , G LL do not involve the boundary, and are therefore determined by the perfect conductance in a Luttinger liquid with coupling parameter z [70], where G 0 = 2z.The heaviside step functions, Θ (± (x − x ′ )), ensure that the responses are causal, and The term G RL corresponds to a non-causal response function, as it measures a right moving current which results from a given left moving current, and therefore vanishes for all x, x ′ , as we verify explicitly in Section IV below.We thus have We can now define the reflection coefficient.The conductance G should be equal to z in the absence of an impurity (E ⋆ J → 0), and vanish for an impenetrable boundary (E ⋆ J → ∞), corresponding to r → 1 and r → −1, respectively.We therefore write The Kubo formula relates the conductance, and hence the elastic phase shift and inelastic decay rate, to a 2point response function.Higher-order response functions yield more refined rates; in this work, we focus on the energy-resolved inelastic decay spectrum, γ (ω ′ |ω), which is the rate of the decay process ω → ω ′ + i ω ′′ i for any choice of ω ′′ i such that ω = ω ′ + i ω ′′ i .The spectrum γ (ω ′ |ω), which may be evaluated in an experiment by injecting a photon ω and measuring the reflected power at frequency ω ′ (see Fig. 1), corresponds to a 3-point response function, as shown in Ref. [31] and discussed in Section V. Measurements of the total inelastic decay rate and spectrum shed light on the fundamental properties of the models.At large enough frequencies, we expect to find Luttinger liquid scaling laws, such as γ (ω) ∼ ω 2z−2 ; the exponents, which could be probed directly in such experiments, indicate whether the boundary operators in ( 4) and ( 8) are relevant (z < 1) or irrelevant (z > 1).In the bsG model, this is of particular importance in the context of the elusive Schmid-Bulgadaev quantum phase transition [47,48], where the relevant and irrelevant cases correspond, respectively, to the superconducting and insulating phases of the boundary Josephson junction.Measuring r (ω) using single photon spectroscopy thus offers a non-invasive approach to probe the dynamics of the boundary junction, avoiding DC transport measurements which were mostly used in previous works [71][72][73][74], and potentially settling a recent debate concerning the presence of the transition [49][50][51][52][53][54].The measurements must be supplemented by theory, notably in the low energy sector, where a perturbative solution does not apply since the boundary cosine in ( 4) is relevant for z < 1. Signatures of the transition could also be found in the spectrum γ (ω ′ |ω) [39].A non-perturbative calculation of the total inelastic decay rate γ (ω), the elastic phase shift δ (ω), and the inelastic decay spectrum γ (ω ′ |ω), for both bsG and Kondo models, is the objective of this work.
An exact solution for the quantities above is possible thanks to the integrability of the models.As explained in the next Section, the Hamiltonians (4) and ( 8) are characterized by an underlying purely elastic dynamics, and may be solved using variations of the Bethe ansatz which apply to field theories.As discussed above, this pure elasticity seems to be at odds with the observed photon splitting.To this end, it should be mentioned that realistic implementations of both models would inevitably suffer from integrability-breaking perturbations.To name a few, stray capacitances, array disorder, and finite lattice spacing, all unavoidable in a cQED environment, would break the integrability of the systems.In addition, the models in Eqs. ( 4) and ( 8) lack a cutoff frequency, which would be present in any physical realization.Therefore, one may be tempted to attribute the inelastic decay of photons observed in cQED experiments to the presence of integrability-breaking terms in the Hamiltonian.However, it is clear in these experiments that inelastic decay persists even when the circuit parameters are pushed towards those of the integrable systems.Particularly, inelastic decay rates observed for a transmon impurity [33] are of the same order of magnitude as those measured for a Cooper pair box [36], where the charging energy of the impurity Josephson junction is pushed to larger and larger values and the Hamiltonian of the system approaches Eq. ( 4).Furthermore, we stress that while experimental systems are not exactly integrable, the use of the integrability formalism is still justified within reasonable assumptions.This is discussed in Appendix B.
Finally, we note that the integrability formalism can only be applied for z < 1 (see explanation below).Fortunately, for z > 1, the RG scale vanishes, and perturbation theory holds for z > 1 at all frequencies (assuming that the boundary Josephson energy is much smaller than the UV cutoff), so that an exact solution is unnecessary.From here on, we restrict ourselves to the regime z < 1.

III. OVERVIEW OF MASSLESS INTEGRABLE QUANTUM FIELD THEORIES WITH A BOUNDARY
In this section, we briefly sketch the main ingredients needed for the calculation of correlation functions in the bsG and Kondo models.These ingredients are the Smatrices, reflection matrices, and form factors.Here we only outline the general structure of the calculation mechanism and introduce some useful notations, relegating explicit expressions to the Appendixes.More details may be found in several papers and reviews [37, 57-59, 75, 76].
It is useful to think of the quadratic bulk Hamiltonian of Eqs. ( 4) and ( 8) as the limiting case of a more general Hamiltonian, For an infinite transmission line (as opposed to the halfinfinite line above), this would be the Hamiltonian of the bulk sine-Gordon model, with a coupling parameter β = √ 8πz, which can be implemented in a cQED environment by connecting the superconducting grains to the ground via Josephson junctions [26].As mentioned above, we focus on the regime z < 1 (β < √ 8π).The nonlinear equations of motion of the generalized bulk Hamiltonian are solved by solitons and antisolitons -wave packets that propagate in the nonlinear medium without dispersion.In the attractive regime, z < 1/2 (β < √ 4π), a pair of a soliton and an antisoliton may form a bound state, called a breather, which also propagates in the bulk without dispersion.The energy and momentum of a soliton (ϵ = +), an antisoliton (ϵ = −), or a breather (ϵ = m, where m is an integer satisfying 1 ≤ m ≤ 1 z − 2), can be parameterized by its rapidity θ as E = M ϵ cosh θ and P = M ϵ sinh θ, where from here on we set the velocity to unity, v = 1.The excitation mass, M ϵ , scales with the bulk Josephson energy, M ϵ ∼ E bulk J 1/(2−2z) [77].We therefore refer to (4) and (8) as massless Hamiltonians; writing θ = ± (A + λ) with A → ∞ such that M + e A /2 → 1, the energy and momentum become those of chiral wave packets, E = µ ϵ e λ and P = ±µ ϵ e λ , where a plus (minus) sign corresponds to a right (left) mover, and µ ϵ = M ϵ /M + (µ ± = 1, µ m = 2 sin (mξ/2)) is the bulk mass ratio.
The presence of an extensive (infinite) number of locally conserved charges in the integrable field theories leads to purely elastic scatterings.Any two solutions to the equations of motion with energies E 1,2 and momenta P 1,2 will maintain their individual energies and momenta when they scatter off one another.The same holds for reflections off the boundary -the momentum and energy of any reflected solution are equal to the incoming ones.This property makes the solitons and breathers a natural basis to diagonalize the Hamiltonians (4) and (8).Quantizing these field theories, a classical field configuration of type ϵ and rapidity θ is a assigned a field excitation, which can be added to a given state by a creation operator Z † ϵ (θ), or removed from a state by an annihilation operator Z ϵ (θ).These elementary excitations are a defining feature of integrable quantum field theories, and form eigenstates of the Hamiltonian, . The dynamics of the excitations is strongly constrained by the extensive number of local conservation laws present in the system, as captured by the Zamolodchikov-Faddeev algebra [9,10]: The S-matrix ) satisfies the Yang-Baxter equation and several additional symmetry relations which are summarized in Appendix C. The relations in Eq. ( 19) are a manifestation of the factorization of nbody interactions to a product of 2-body interactions in the integrable system.
The eigenstates |θ n , . . ., θ 1 ⟩ ϵn...ϵ1 form a complete set of states.In the massless limit and in the presence of a boundary, the complete set is [59] where the states take into account both right and left movers: The superscript a of a rapidity indicates the momentum carried by the excitation -it is ς a µ ϵ e λ , where ς r = 1 and ς l = −1.Here R ϵ ′ ϵ (λ) is the boundary reflection matrix, which is non-zero only if µ ϵ = µ ϵ ′ ; therefore, a soliton can become an antisoliton (or vice-versa) upon reflection of the boundary, but a reflected breather m remains a breather of the same type.As shown in [58], R ϵ ′ ϵ (λ) depends on the difference λ − λ B , where λ B is the rapidity associated with the energy scale of the boundary, ), and given explicitly by [78] where Γ (x) is the gamma function [79], and ε is either E J or J for the bsG or Kondo Hamiltonians, respectively (recall that Λ is the UV cutoff frequency).Note that the reflection matrices are thus defined only for z < 1, where T B is finite; hence, this formalism cannot be used to treat devices with z > 1.The reflection matrices for the bsG and Kondo models are given in Appendix D. In Eq. ( 21), Sϵ1ϵ2 (θ) denotes the S-matrix for the exchange of a right mover and a left mover; note that it is diagonal.
The complete set of states can be inserted between any two operators in a correlation function.For example, consider a 2-point correlator: In the following, if the coordinates of the operator are omitted, it should be understood that the matrix element is evaluated at the origin.Before we proceed, let us introduce a few additional notations that will simplify the equations to follow., where ± = ∓ and m = m.The complete set of states in Eqs.(20), ( 21) is then compactly written as where and the product of reflection matrices is abbreviated as The form factors are appropriately abbreviated as In many of the following calculations, the rapidities will be shifted according to their bulk mass ratios µ ϵ = M ϵ /M + : λ k → λ k − log µ ϵ k .This will be denoted by a hat: and similarly for scattering and reflection matrices.We also define the total energy of the shifted set ⃗ λ, ν λ ≡ ) and, using the abbreviated notation, The boundary is depicted by Note that the order of the excitations is reversed upon reflection off the boundary.In the abbreviated notation, The form factors in the massless sine-Gordon model are given in Appendix E. The form factors, scattering matrices, and reflection matrices are all of the ingredients we need to calculate our desired response functions.

IV. THE TOTAL INELASTIC DECAY RATE AND ELASTIC PHASE SHIFT
We are now equipped with all of the necessary tools to calculate response functions.In this section, we calculate the reflection coefficient r (ω), which defines the total inelastic decay rate γ (ω) and the elastic phase shift δ (ω) (Eq.( 11)), using Eqs.( 17) and (14).

A. The AC conductance
The conductance G LR (x, x ′ ; ω) has been calculated in Ref. [59] for the boundary sinh-Gordon and bsG models.The conductance of the bsG model with z = 1/3 is of particular interest, as it describes the tunneling of fractionally-charged excitations in a fractional quantum Hall sample at filling ν = 1/3.Closely-related expressions were also obtained for the Kondo model.Here we retrace the steps of Ref. [59] and calculate G LR (x, x ′ ; ω) for the bsG and Kondo models with general coupling parameters.We emphasize the key steps in the calculation of the 2-point response function, laying the groundwork for the derivation of the 3-point response function in the next Section.We also show that the non-causal response function G RL (x, x ′ ; ω), measuring the response of a right moving current to a left moving current perturbation, is manifestly zero in all models, for any coupling parameter and impurity strength, as it should be.
The Kubo formula relates the conductance G LR (x, x ′ ; ω) to a 2-point response function of the current operators.Each of the two correlators in Eq. ( 14) can be evaluated via a form factors expansion, by inserting the complete set of states in Eq. ( 25).As we will now show, it is necessary to combine the two correlators to obtain a single unified expansion for G LR (x, x ′ ; ω).Let us start with ⟨L (x, t) R (x ′ , 0)⟩; using Eqs.( 25), ( 24), (E16), and (E19), we find Shifting the rapidities according to the bulk mass ratios of their corresponding excitations, We now introduce an auxiliary rapidity, κ: Shifting the rapidities by κ, λ k → λ k + κ, we obtain where we used Lorentz invariance, Eq. (E3), and the notation ⃗ λ + κ = {λ 1 + κ, . . ., λ n λ + κ}.For the second correlator of the commutator, ⟨R (x ′ , 0) L (x, t)⟩, follow-ing the same steps, we find Now, we shift κ → κ − iπ in Eq. (37).Using Next, plug Eqs. ( 36) and (38) to Eq. ( 14).Integrating over time yields The integration contour of κ can be closed by adding the edges at κ = σ + iχ, with σ → ±∞ and 0 ≤ χ ≤ π, where the integrand decays rapidly (recall that x, x ′ < 0).In the presence of bound states, z < 1/2, there could be poles of the reflection matrices enclosed by the contour; however, the residue of the integrand at κ = κ 0 gives rise to the exponential e e Re{κ 0 } sin(Im{κ0})(x+x ′ ) , where sin (Im {κ 0 }) > 0, hence this term vanishes in the limit of a half-infinite line for large enough |x + x ′ |.Therefore, the only pole contributing to the integral is at κ = log ω + iη: Before we proceed to the calculation of the reflection coefficient, it is instructive to consider the non-causal conductance G RL (x, x ′ ; ω), and show that it vanishes.An identical treatment would lead to Note the three crucial differences compared to Eq. ( 39): the complex conjugation of the reflection matrices, the sign in the space-dependent exponential, and the integral dκ, which runs along the line Imκ = −iπ.This time, the pole at κ = ω+iη is not enclosed by the contour, leading to G RL (x, x ′ ; ω) = 0.

B. The reflection coefficient and the origin of inelastic decay
We now return to Eq. ( 40) and extract the reflection coefficient using Eq. ( 17).Reshifting the rapidities as k=1 e λ k is the sum of energies of the excitations in the set λ.The physical interpretation of the expression above is clear -the reflection coefficient of a photon with frequency ω is given by the sum over all excitations with energies summing up to ω.Note that this result is general, and applies to any integrable boundary model -the choice of model specifies the reflection matrices and form factors.The only energy scale, T B , is encoded within the reflection matrices.While the sum in (42) runs over all possible number of excitations n λ , the terms decay very rapidly with n λ , and in practice it is enough to calculate only a few terms to obtain r (ω) with excellent accuracy.The accuracy of the expansion can be evaluated using the high frequency behavior of r (ω); using )), we find r (ω ≫ T B ) → ϵ λ r 0 ϵ λ = 1 in all models, where Then, calculating the truncated sum ′ ϵ λ r 0 ϵ λ provides an estimation for the accuracy of the result, and sets an upper bound on the contribution of the remaining terms (since r 0 ϵ λ > 0 for all ϵ λ ).At low frequencies, we have r for most values of z.The leading terms for the bsG and Kondo models are as follows.First, a soliton-antisoliton pair contributes This is the leading term for z ≥ 1/2.In the attractive regime, z < 1/2, the dominant contribution comes from a single breather: where and In the following, we evaluate the reflection coefficient at integer p = 1/z as The evaluation of r +−1 is significantly more complicated for non-integer 1/z.However, Fig. 4 shows that 12 becomes larger than that of r R +−1 ), and also for z > 1/2, as long as z is not too large (where r +−+− is non-negligible).The inelastic decay rate and elastic phase shift are then readily extracted following Eq.(11).
Eq. ( 42) reveals the origin of the inelastic scattering.The inelastic decay rate of a photon ω is We thus identify the coherent sum in Eq. ( 42) as the source of the photon decay.A plane wave mode at frequency ω impinging on the boundary can be formally written as a sum of eigenstates, | ⃗ λ⟩ r ϵ λ , with appropriate weights and phases.Note that this is a nonlinear relation -plane waves cannot be expressed as a sum of individual solitons and breathers, and must be spanned using all eigenstates with any number of excitations.While each excitation in each of the eigenstates is reflected elastically off the boundary, it picks up a phase, determined by the reflection matrix, that depends on its type and energy.These relative phases between the eigenstates alter the specific weights of the decomposition of the incoming photon.Therefore, the reflected excitations no longer form a single photon at frequency ω, but rather a set of photons with frequencies ω i , such that i ω i = ω.
The nonlinear relation between the photons and elementary excitations is generally implicit, and is hidden within the form factors of the theory, which are the matrix elements of the derivatives of the bosonic field ϕ in the basis of solitons and breathers.An explicit relation can be found at the free-fermion point, z = 1/2, using refermionization [80], as discussed in Appendix G. Introducing a fermionic field ψ (x) ∼ e iϕ(x) , we find a oneto-one correspondence between ψ, ψ † and solitons and antisolitons by expanding the fermionic field as k satisfy the Zamolodchikov-Faddeev algebra (Eq.( 19)) at z = 1/2, where the S-matrix becomes trivial, We therefore identify ψ k>0 with a soliton with energy k, and ψ † −k<0 with an antisoliton with energy k, establishing an explicit nonlinear relation between the solitons and the bosonic field ϕ.The calculation of r (ω) simplifies considerably at the free-fermion point, as the sole contribution to the reflection coefficient is Eq. ( 44), and the particularly compact expressions for the reflection matrices (Eqs.(D3)-(D6)) and form factors (Eq.(E32)) allow for closed analytical expressions: The same expressions are derived in Appendix G using refermionization.

C. Results
The rates γ (ω) and δ (ω) in both models, for several values of z, are displayed in Fig. 2. The inelastic rate in both models follows a Luttinger liquid power law at high frequencies, γ (ω ≫ T B ) ∼ ω 2z−2 , which may be obtained from the Hamiltonians (3) and ( 5) by means of perturbation theory [31,35,36].Note that Eq. ( 22) restores the dependence of the rates on the UV cutoff Λ (as long as ω ≪ Λ), which is present in a perturbative treatment and should be considered for quantitative comparisons with experimental measurements, not only in the high frequency regime, but also in the low frequency regime, ω ≪ E ⋆ J .Indeed, perturbation theory is invalid for the Hamiltonians (3) and (5) for z < 1 and ω ≪ E ⋆ J .The low frequency power laws for the rates are model dependent, and could be predicted from the expansions near the strong coupling fixed points.In the Kondo model, the dominant contribution stems from a quartic density term, ρ 4 , leading to γ K (ω ≪ T B ) ∼ ω 6 for all z < 1 [31].In bsG, the leading expansion terms are the quartic phase term ϕ 4 , giving rise to the same ω 6 power law as the quartic density term, and the dual cosine cos (πρ), which generates instantons between the minima of the cosine potential, and leads to a ω 2/z−2 behavior [40]; hence, γ bsG (ω ≪ T B ; z ≥ 1/4) ∼ ω 2/z−2 and γ bsG (ω ≪ T B ; z < 1/4) ∼ ω 6 .
Figure 2. Total inelastic rate γ (ω) and phase shift δ (ω) for the bsG and Kondo models and several values of z.The power laws at low and high frequencies are denoted on the plots.In the phase shift panels, we plot both δ (ω) and π/2 − δ (ω) or π − δ (ω) (for the bsG and Kondo models, respectively).We use Eq. ( 48) to evaluate r (ω) at integer p = 1/z, and The asymptotic behavior of the phase shift δ (ω) may be obtained from Eq. ( 42), rewritten here as As discussed in Ref. [59], one may expand the reflection matrices as a power series in ω/T B or T B /ω, for low or high frequencies, respectively.In the bsG model, the soliton matrices R ±/∓ ± (λ) expand as double power series in e λ and e (1/z−1)λ , and the soliton matrices R m m (λ) expand as power series in e λ .We find δ bsG (ω ≫ T B ; z ≥ 1/2) ∼ ω 2z−2 and δ bsG (ω ≫ T B ; z < 1/2) ∼ ω −1 at high frequencies (in agreement with perturbation theory), and π/2 − δ bsG (ω ≪ T B ; z ≥ 2/3) ∼ ω 2/z−2 and π/2 − δ bsG (ω ≪ T B ; z < 2/3) ∼ ω at low frequencies, which are again anticipated by the expansion near the strong coupling fixed point.The reflection matrices of the Kondo model are considerably simpler and all expand as power series in e λ , leading to We note that all of the scaling laws of the scattering rates γ bsG (ω) and δ bsG (ω) in the bsG model, at frequencies both above and below the RG scale, are in agreement with a recent theoretical study combining perturbation theory at ω ≫ T B and strong coupling expansions at ω ≪ T B [42]; our exact calculation provides the full crossover between the two regimes.
To conclude this section, we relate the limiting cases of the inelastic and elastic rates to the Schmid-Bulgadaev transition.Approaching the transition, z → 1, we find that both low and high frequency power laws, 2/z −2 and 2z − 2, respectively, tend towards 0. In other words, the rate γ bsG (ω) and phase shift δ bsG (ω) become frequencyindependent at the transition point, z = 1.Perturbation theory at z > 1 shows that γ bsG (ω) ∼ ω 2z−2 for both low and high frequencies (since the boundary cosine term of Eq. ( 3) becomes irrelevant at z > 1), and δ bsG (ω; 1 < z ≤ 3/2) ∼ ω 2z−2 and δ bsG (ω; z > 3/2) ∼ ω, also for all frequencies [36].Crucially, the phase shift at ω = 0 jumps from π/2 to 0 across the transition, and changes sign at finite frequencies.The sign change of δ bsG (ω) and the different trends of γ bsG (ω) on the two sides of the transition, both at finite frequencies, were recently observed using single photon spectroscopy, providing evidence for the long-sought-after quantum phase transition [36].

A. The spectrum as a 3-point response function
The inelastic decay spectrum, γ (ω ′ |ω), measures the production rate of photons at frequency ω ′ due to a splitting of a photon with frequency ω > ω ′ .It is related to the total inelastic decay rate, γ (ω), by means of an energy-conservation sum rule: It was shown in Ref. [31] that, for the generic impurity setup in Eq. ( 9), γ (ω ′ |ω) can be found by calculating a 3-point response function: which is defined in the time domain as where x in < 0 is an arbitrary point in the half-infinite lead, and n q ′ = b † q ′ b q ′ is the occupation number of the mode q ′ , with frequency ω ′ = vq ′ = q ′ .The Fourier transform is Let us sketch the main steps leading to this expression.The goal is to calculate the time-averaged rate of change of n q ′ in response to an incoming photon with frequency ω, injected at x in < 0, impinging at the boundary.Note that we are looking for a second order response, since the change in the photonic number n q ′ is proportional to the flux of incoming energy.The photon at ω is injected by applying an AC voltage, H AC = V (t) e ηt ρ (x in ), where V (t) = 2V 0 cos (ωt) and η → 0 + .The time dependence of H AC makes the Keldysh formalism natural for our purpose [81,82]; the second order Kubo formula for n q ′ reads ⟨n (57) The notation c, q corresponds to the "classical" and "quantum" fields in the Keldysh formalism, that is, the sum and difference, respectively, of the fields on the forward and backward time contours.We are concerned with the so-called "fully-retarded" multipoint correlator G cqq ω ′ [83], which measures the causal response of n q ′ to the perturbing density ρ.To find the time-averaged photon production rate, we take the derivative of Eq. ( 57) and discard the oscillating exponentials: Then, in order to obtain γ (ω ′ |ω), we multiply the above by the density of modes, 1/∆ = ℓ/π, and divide by the average power of the source that propagates towards the impurity, ω V 2 0 /π.We thus arrive at Eq. ( 54).Using time-translation invariance and applying simple algebraic manipulations, we may write the 3-point response function as Note the η → 0 + prefactor in Eq. ( 54), which implies that we need to look for contributions to G cqq ω ′ that are singular in η.
In order to calculate G cqq ω ′ using form factors, we need to express n q ′ , ρ in terms of the current operators R, L. From Eqs. (10), (13), and the quantization q = n∆, we may write b where ω = vq = q, ς R = 1, and ς L = −1.The operators b R q and b L q annihilate right and left moving photons with frequency q, respectively.The density ρ and occupation number n q ′ may then be decomposed into their chiral parts.The density can be written as ρ = ρ R + ρ L , where whereas the occupation number may be expressed as The decomposition of G cqq ω ′ into its chiral parts is then with Fortunately, there is no need to calculate all 16 terms G cqq ω ′ ;ABCD .First, in an experimental setup, the photon at ω is injected such that it propagates towards the boundary; we should therefore only consider terms with ρ C , ρ D = ρ R .This leaves us with the calculation of 4 terms, G cqq ω ′ ;ABRR .Furthermore, since the produced photon propagates to the left, away from the impurity, we expect only G cqq ω ′ ;LLRR to contribute, measuring the response of the left moving occupation number to the right moving source.In the following, we calculate G cqq ω ′ ;LLRR , and show that the other 3 terms indeed vanish.

B. Calculating the 3-point response function using form factors
The calculation of a 3-point response function using form factors is considerably more involved than that of a 2-point function.The 3-point response function G cqq ω ′ ;LLRR comprises 4 correlators, each involving 4 current operators.We therefore must insert 3 complete sets of states in each correlator, giving rise to mixed matrix elements of the form ⟨ ⃗ ϑ|A| ⃗ θ⟩ ϵ ϑ a a ϵ θ .These mixed elements, evaluated using the crossing relations in Eq. (E11), lead to a series of terms with a different structure in each of the 4 correlators.
In order to make sense of these complicated expressions, let us recall the key step in the derivation of the 2point function G LR (t).It comprises only two correlators, ⟨L (t) R (0)⟩ and ⟨R (0) L (t)⟩, which were both evaluated by inserting a single complete set of states.The result, however, was not obtained by considering each correlator on its own, but rather by combining the two form factor expansions; we had to take the difference of these two correlators, allowing us to close the integration contour in Eq. (39), in order to arrive at the anticipated delta function in the result, Eq. ( 40).Similarly, it was necessary to combine the two correlators of G RL to show that the non-casual conductance vanishes.It is therefore crucial to find some convenient way to combine and unify the 4 correlators of G cqq ω ′ ;LLRR ; given a term in the form factor expansion of one of the correlators, we must find a way to identify its 3 counterparts in the other correlators.In the following, we accomplish this and show how to identify a quartet of terms from the expansions of the four correlators that need to be summed up together, by labeling the excitations according to their "origins" and "destinations" in each correlator, leading to a unified general expression for G cqq ω ′ ;LLRR .We begin by denoting the 4 terms of the double commutator as so that Consider III above and insert complete sets of states between its operators: Using Eq. ( 24), we get complex exponentials from the coordinates of the current operators.We assume a half-infinite line and extend the lower integration limits of x 1,2 to −ℓ → −∞, introducing an infinitesimal parameter η to assure convergence: Integrating over x 1,2 and t ′ and shifting the rapidities according to the bulk mass ratios of their excitations, as was done in the calculation of the reflection coefficient, then yields where a hat over the operators R, L indicates that the rapidities in the matrix elements are λ − log µ ϵ .Now, consider the mixed matrix elements, which need to be evaluated using the crossing relations in Eq. (E11).The idea is to partition the sets λ 1,2,3 into smaller subsets, α ij (i, j = 1, 2), β, and γ, that label the excitations according to the operators they are connected to: the excitations in the set to each other, and the set γ connects the ρ R operators.We find The derivation of Eq. ( 70) using the crossing relations is detailed in Appendix F.Here we shorten the notations, Counting the number of possibilities to partition the sets λ i to α ij , β, γ, we see that the integration measure becomes that is, the factorials and 2π factors in the denominators translate naturally from the λ sets to the smaller subsets.The other correlators, I, II, IV, may be treated similarly, with appropriate labeling of the excitations.The 4 correlators, written explicitly in Eqs.(F4)-(F7), are conveniently represented by the graphical notation: Labeling the excitations offers a natural way to sum up the correlators, following a procedure similar to that of the 2-point function.Consider I, III; we introduce two auxiliary rapidities, κ 1,2 , such that e κi = ν i1 + ν i2 , and shift where C λ is a sign factor stemming from charge conjugation, defined in Eq. (E9).Shifting Note that we also rewrote the ω ′ denominators using the delta functions.In order to close the κ 1 contour, we add the edges at κ 1 = σ + iχ with σ → ±∞ and −π ≤ χ ≤ 0 -this is allowed, since the exponential e i(e κ 1 −e κ 2 )xin assures that the integrand decays fast enough at σ → ∞, whereas the form factors vanish exponentially-fast as σ → −∞.Now, the pole at κ 1 = κ 2 − 2iη is enclosed by the contour; setting κ 1 = κ 2 , the ω ′ denominators become that is, the contribution of this pole is singular in η.In fact, the pole at κ 1 = κ 2 − 2iη is the sole singular contribution to the contour integral.The annihilation poles of the form factors, occuring at θ i − θ j = iπ for some i > j (when the same excitation appears on both bra and ket of a matrix element, see Eq. (E4)), are just above or below the upper and lower boundaries of the contour, and while in the attractive regime (z < 1/2) the form factors and reflection matrices have bound state poles that are enclosed by the contour, their residues are not singular in η.
We are now in position to show why G cqq ω ′ ;LLRR is the only non-vanishing contribution to the spectrum.First, G cqq ω ′ ;RRRR is a background term that does not involve the boundary and thus cannot account for inelastic scattering.Next, consider G cqq ω ′ ;RLRR , that is, the response of n RL q ′ to a right-moving photon.Replacing in Eq. ( 67), the sign in the x 1 exponential in Eq. ( 68) is flipped: The key difference with respect to Eq. ( 67) is the same relative sign of the iη terms in the denominators.Therefore, G cqq ω ′ ;RLRR (and, similarly, G cqq ω ′ ;LRRR ) is non-singular in η, and its contribution vanishes in the limit η → 0. It is also reassuring to verify that G cqq ω ′ ;RRLL , the non-causal counterpart to G cqq ω ′ ;LLRR , vanishes identically; this is shown in Appendix F.
The treatment of II, IV is identical to the above.Plugging everything to Eq. ( 66) and back to Eq. ( 54), we find where the products of form factors and reflection matrices are (denoting Note that we subtract the product of Kronecker deltas in R.This corresponds to subtracting the background term G cqq ω ′ ;LLLL , which does not involve the boundary and therefore does not contribute to inelastic scattering, and assures that γ (ω ′ |ω) vanishes as ω/T B → ∞ (following Eq. (D8)) or ω/T B → 0 (following R where P denotes the principal value.To get rid of the awkward principal value terms, consider another equivalent way to expand the correlators I, . . ., IV: Note that the diagrams III and IV above are the mirror images of the ones in Eqs. ( 74) and (75), where bra states become ket states and vice-versa -this corresponds to complex conjugation of both form factors (see Eq. (E14)) and reflection matrices.This time, we pair I with IV and II with III.The same steps lead to The expressions in Eqs. ( 81), ( 89) are equal; taking their average eliminates the term proportional to Im { F R}, leaving us with the simple delta function δ (ω − Ω).Finally, we arrive at This is a general form factor expansion for the inelastic spectrum.
Before we proceed, let us consider the disconnected case, ϵ 12 = ϵ 21 = ϵ β = ϵ γ = {}, for which we find (91) That is, the disconnected terms do not contribute to the spectrum at ω ′ < ω, and reproduce the total inelastic decay rate.That is similar to the calculation in the Keldysh formalism, where the total rate is also obtained from the disconnected diagrams (see Appendix G and Ref. [31]).

C. Leading diagrams
Eq.
(90) offers a non-perturbative and rapidlyconvergent diagrammatic approach for calculating the spectrum.Each term in Eq. ( 90) may be represented by a diagram, from which one may read off the corresponding form factors and reflection matrices, as illustrated below.The physical intuition behind this approach is clear; one needs to sum over all processes with excitations of the integrable theory connecting the bosonic operators, with energy conservation imposed by the delta functions.
Comparing to refermionization at the free fermion point z = 1/2, discussed in Appendix G, we draw an analogy between our approach and Wick's theorem, as explained below.
All that is left to be done is to draw the leading diagrams and evaluate their contributions.As in any form factor expansion, such as that of the reflection coefficient in Eq. ( 42), we expect the contributions to decay rapidly with the number of excitations involved, allowing us to obtain the spectrum with good accuracy using only a few terms.Our figure of merit is the sum rule in Eq. ( 53), which should hold for all ω; hence, it is essential that the asymptotic power laws of γ (ω) at low and high frequencies will be recovered by the sum In the following, we focus on z = 1/3 and z = 1/2, list the leading terms in Eq. ( 90), and draw their corresponding diagrams.The details behind the evaluation of these terms, as well as additional subleading terms, are given in Appendix F. In particular, one must be careful to take into account all terms with comparable contributions, which could be misleading in the presence of mixed matrix elements, ⟨ ⃗ ϑ|A| ⃗ θ⟩ ϵ ϑ a a ϵ θ , due to the annihilation poles of the form factors (Eq.(E4)).Appendix F presents a consistent method to identify and evaluate all such terms.
First, consider the diagram s′′ .
(92) Its contribution to the spectrum reads with e λ1 = ω −Ω, e λ2 = Ω, e λ3 = ω −ω ′ −Ω, e λ4 = ω ′ +Ω, as specified on the diagram lines.Note that we take the real parts of F and R separately; this follows from Eq. (F18) (see Appendix F for details).At the freefermion point, z = 1/2, all form factors other than f R +− vanish (see Appendix E), and we find that this is the only contribution to the spectrum, γ (ω 1 (ω ′ |ω).As shown in Fig. 3, the sum rule is indeed perfectly obeyed in that case, for both bsG and Kondo models.Furthermore, we show in Appendix G that the same results for either model can be obtained by means of refermionization.There is a clear connection between Eq. ( 92), where each ρ R leg is connected to both b L † q ′ and b L q ′ , and the contractions in Eq. (G16) that lead to Eq. (G22).This agreement not only highlights the one-to-one correspondence between the solitons and fermions at z = 1/2, but also draws an analogy between our diagrammatic approach and Wick's theorem in free theories, and serves as an essential sanity check for the general expression in Eq. (90).
Varying z in the vicinity of z = 1/2, we expect γ (1) 1 (ω ′ |ω) to dominate; indeed, for z ∼ 1/2, the sum rule is obeyed with good precision.The precision deteriorates as we drift away from the free-fermion point, where the contributions of multi-soliton (z > 1/2) or breather (z < 1/2) form factors become important, and considerably expand the available phase space beyond the diagram in Eq. ( 92).However, panel (b) of Fig. 3 shows that the z-dependent power laws of the sums ω 0 ω ′ ω γ (ω ′ |ω) dω ′ are equal to those of γ (ω); a power law mismatch would be reflected by a sharp increase or decrease of the ratio at low or high frequencies.
Moving on, we concentrate on z = 1/3, where form factors involving the m = 1 breather should play a major role.We find the leading term to be where e λ1 = ω, e λ2 = ω − ω ′ − Ω, e λ3 = Ω, e λ4 = ω − Ω.Again, we take the real parts of both F and R, as detailed in Eq. (F14).In γ In addition to the two diagrams above, we consider 3 additional diagrams in the evaluation of the z = 1/3 spectrum, listed in Appendix F. We find that, for z = 1/3, the resulting spectrum obeys the sum rule with good accuracy for a wide range of frequencies ω, as shown in Fig. 3. Furthermore, the contributions decay rapidly, as demonstrated in Fig. 5 in Appendix F. It is important to note that plugging each of the terms γ i (ω ′ |ω) into the sum rule recovers the correct power laws of γ (ω), for both ω ≪ T B and ω ≫ T B .

D. Results
The spectra for both models at z = 1/3, 1/2, as well as the validation of Eq. ( 53) for several values of z, are displayed in Fig. 3.We have seen in Section IV that the bsG and Kondo models exhibit the same Luttinger liquid power laws at high frequencies, ω ≫ T B , but differ below the RG scale, ω ≪ T B .The spectrum, determined from a higher order response function, gives us more refined information.Indeed, the difference between the models is emphasized by the spectrum; while the Kondo spectrum is suppressed at low produced frequencies ω ′ ≪ T B , the bsG spectrum diverges, γ bsG (ω ′ ≪ T B |ω) ∼ ω ′−1 , regardless of the incoming frequency ω (note that the total rate, given by the sum rule in Eq. ( 53), is still finite).This asymptotic behavior results from the term γ (1) 1 (ω ′ |ω) in Eq. ( 93) -one may show that the product of reflection matrices is different than 1 upon setting ω ′ = 0, and the ω ′−1 behavior results from the prefactor of the integral.In other words, the production of low frequency photons is favored in a splitting process in the bsG model.The proliferation of low frequency produced photons results from the inductive coupling of the halfinfinite line to the impurity Josephson junction in Eq. ( 1), with HI = −E J cos φ0 .The sole source of photon splitting is the boundary cosine term; fluctuations of the flux φ0 in this nonlinear potential give rise to the inelastic decay.The inductive coupling to the line translates fluctuations of φ0 into approximately uniform shifts of φn>0 , which correspond to low frequency modes.This behavior is opposed to the Kondo spectrum, which vanishes at ω ′ = 0; expanding the product of reflection matrices in Eq. ( 93) at ω ′ ≪ T B , we find γ K (ω ′ ≪ T B |ω) ∼ ω ′ .This asymptotic behavior turns out to hold for all contributions to the spectrum and hence to all z, in agreement with Ref. [31].
It is also interesting to consider decay processes with ω ′ ≲ ω, where the total energy of the remaining photons is small, ω − ω ′ ≪ T B .One may then apply strong coupling fixed point expansions to extract the power law dependence of the spectrum as a function of ω − ω ′ .Analytical expressions for the spectrum, and therefore for the ω − ω ′ dependence, may be calculated at the free-fermion point, z = 1/2 (see Eq. (G25), which was obtained using refermionization and is identical to the form factors result).For the Kondo spectrum, one finds 3 , in agreement with the contribution of the strong coupling fixed point expansion of Ref. [31].In the bsG model, one finds γ bsG (ω ′ ≲ ω|ω) ∼ ω − ω ′ at z = 1/2; interestingly, this power law dependence does not stem from neither the quartic term of the boundary cosine operator, nor from the dual cosine term, which both induce a (ω − ω ′ ) 3 power law [40].For z < 1/2 it appears difficult to numerically evaluate the involved integral to sufficient accuracy to determine the behavior with enough certainty, since strong cancellation occurs between different diagrams.We stress that the spectrum at z = 1/2 was obtained, for both models, in two independent methods (form factors and refermionization), which yield identical results.Furthermore, as evident from Ref. 3, the asymptotic power laws of the total decay rate are recovered by the sum rule in Eq. ( 53) for all values of z.It could be interesting to look further into this issue in a future work.

VI. CONCLUSIONS
In this work, we have shown how inelastic decay of microwave photons, measured in cQED experiments implementing integrable systems, can emerge from the purely elastic scattering of the excitations from the integrability picture.Using the framework of form factors, we identified the origin of the photon splitting as the nonlinear relation between the microwave photons and the elementary excitations of the bsG and Kondo models.The form factors, encoding this nonlinear relation, allowed us to obtain exact results, going beyond previous perturbative calculations [35,36,[39][40][41][42]. Crucially, our results hold at low energies, even if the boundary impurity terms of the Hamiltonians are relevant, rendering perturbation theory invalid, as well as when strong coupling expansion fails.The low energy results for the total inelastic decay rate and the elastic phase shift distinguish between the bsG and Kondo models, which both exhibit Luttinger liquid power laws above the RG scale.This distinction is emphasized by the energy-resolved inelastic decay spectrum, where γ bsG (ω ′ ≪ T B |ω) ∼ ω ′−1 and γ K (ω ′ ≪ T B |ω) ∼ ω ′ for all z; note that this result could not be obtained using perturbation theory, even for ω ≫ T B .As discussed in Section IV, such exact low energy expressions are particularly useful in the context of the Schmid-Bulgadaev quantum phase transition, and shed single-photon light on this intriguing phenomenon.
In the calculation of the energy-resolved inelastic decay spectrum, we have devised a general method to calculate a 3-point response function using form factors.While previous works have dealt with the calculation of multipoint correlation functions in integrable quantum field theories using form factors [61][62][63], they focused on theories with a single excitation type and only considered specific contributions.Our method, in contrast, provides a general expression for all orders, and for any excitation content of the theory.The physical intuition behind this diagrammatic approach is clear -one has to sum over all combinations of excitations connecting the bosonic operators, imposing appropriate energy conservation.As in any form factor expansion, it is sufficient to consider a few terms to obtain a result with good precision, numerically evaluating only single or double integrals.The generalization of our method to 4-point response functions or higher is straight-forward.
Looking forward, there are several possible extensions to this work which could improve quantitative comparisons with experimental measurements.The most pressing issue is the incorporation of finite temperatures into our framework; indeed, realistic temperatures in cQED experiments, T ∼ 50 mK ∼ 1 GHz, are usually larger than the RG scale defined by the impurity.In those cases, a perturbative approach is valid at any frequency, and diagrammatic techniques have been shown to provide results that quantitatively agree with experiments [35,36].Yet, the rapid rate of technological advancements in the field of cQED indicates that T ≲ E ⋆ J could soon become possible, thus raising interest in exact finite temperature results at all frequencies.We derive such exact results for z = 1/2 using refermionization in Appendix G (see also Ref. [36] for explicit expressions); other values of z have to be treated within the framework of form factors.In fact, the calculation of finite temperature correlators using form factors [84][85][86][87][88][89][90] involves mixed matrix elements of the form ⟨ ⃗ ϑ|A| ⃗ θ⟩ ϵ ϑ a a ϵ θ , much like those appearing in the 3-point response function considered in this work; hence, the methods applied here could also be useful for evaluating finite temperature response functions.Two other aspects which should be addressed are the finite volume of the system [91,92] and the presence of integrability-breaking terms [93], both of which are inevitable in experimental setups.Furthermore, inelastic decay in cQED setups should not necessarily emerge from reflection off impurities, and can also occur in nonlinear bulk models.Our framework could then be extended to the massive bulk sine-Gordon model, amenable to realization in superconducting circuits [26], and applied to investigate bulk effects, such as the superconductor to insulator phase transition in an array of Josephson junctions [28,[94][95][96], through the lens of elastic and inelastic scattering.
Finally, our results should pave a path for tackling other systems, beyond the scope of cQED experiments.Our framework for calculating form factor expansions of multipoint response functions could be used in other contexts of integrable field theories; a particular example is the calculation of 4-point functions in tunneling experiments between fractional quantum Hall leads at filling ν = 1/3, whose low energy behavior is captured by the bsG model with z = 1/3.Another likely application is in the context of one-dimensional cold atom systems [97], which are often integrable.In particular, we could use our developed formalism to evaluate multipoint response functions that are measured in post-quench evolutions and indicate non-gaussian correlations which are ubiquitous in integrable systems [18,98,99].Multipoint functions could also be used to investigate the onset of chaos due to weak integrability-breaking terms, which, when treated perturbatively within the form factors formalism, necessitate the use of higher-order correlation functions the order of the desired correlation functions [100][101][102].Our method could also be applied for the calculation of multipoint correlators in high energy contexts, such as the N = 4 supersymmetric Yang-Mills theory in 3 + 1 dimensions [103], which is suspected to be integrable.
with the impurity Hamiltonian where Φ ext is an external magnetic flux, which we take from here on as half flux quantum, Φ ext = πℏ/ (2e) = π/2.It is useful to rewrite HI as We proceed by diagonalizing the array Hamiltonian Ha .Hamilton's equations read where ω ∥ = α/ C f L ∥ .We look for an oscillatory solution, φn ∼ k φk e −iω k t .In order to decouple the φ0 term from the n > 0 equations, we define new variables, φn = φn − ω 2 ∥ φ0 /ω 2 k , leading to In the following, we set the array spacing to unity, a = 1.The equations are solved by φn ∼ sin (kn − δ k ).The bulk equations, n > 0, yield the dispersion relation, ω k = 2v sin (k/2) ≈ vk, where the velocity v = 1/ C g L line is assumed to be much larger than any other energy scale.
The n = 0 equation yields the phase shift, L line /C g is the inverse RC time of the fluxonium and the transmission line.The approximation above holds provided that α 2 C g /C f ≪ 1 and L ∥ ≫ L line , so that 1/ L ∥ C g is significantly larger than all energy scales other than v.We impose open boundary conditions at n = N , leading through the Hamilton equation for n = N to sin (kN − δ k ) = sin (k (N + 1) − δ k ), which yields a quantization condition kN − δ k = πm + π/2 with m = 0, 1, . . ., N , and therefore a mode spacing ∆ ≈ πv/N .In the following, we also need the capacitance matrix [C] n,m , obtained by inverting the capacitance energy matrix of Ha : Neglecting 1/N corrections, we find that the mode capacitances (that is, the eigenmode expectation values of the capacitance matrix) are given by The line Hamiltonian may now be quantized by introducing creation and annihilation operators.The diagonalized phase and charge operators read and the diagonalized array Hamiltonian is given by Using the quantization condition, we find allowing us to express the coupling term H c in terms of the eigenmodes: (A11) The square root in the denominator above imposes a high-frequency cutoff.To make contact with standard bosonization expressions, we replace it by an exponential cutoff, e −ω k /Λ , where Λ ∼ min ω ∥ , Γ f .Finally, we are in position to derive the spin-boson Hamiltonian.The fluxonium may be approximated as a symmetric double well potential for the flux ϕ f (achieved by tuning the external magnetic flux to half flux quantum [65]).We label the two lowest eigenstates of the fluxonium Hamiltonian as the eigenstates of the S x operator, so that S z eigenstates correspond to wavefunctions with well-defined phase, localized near either of the minima of the double well.The fluxonium Hamiltonian then reads H f = −JS x , where J is the tunneling matrix element between the two wells, which can be calculated by WKB or instanton methods [104].The fluxonium charge Qf couples the two-level system to the array via a S z term [66], where is the expectation value of the fluxonium charge operator in the eigenfunctions of either of the symmetric wells.The full Hamiltonian reads This k-space version of the spin-boson Hamiltonian can now be easily mapped to the continuum real-space version of Eq. ( 6).Note that the coupling coefficient is given by √ 2z ′ πv, with z ′ = z × (2q f /α) 2 ; that is, the Luttinger parameter of the Kondo model is proportional, but not equivalent, to the normalized impedance z.
where m = 0, 1, 2, . . .and δ 0 is a phase shift set by the boundary conditions at the far end of the line, away from the impurity.Our analysis should hold as long as we are concerned with modes at frequencies ω ≫ ∆; note that in realistic setups [33][34][35][36], the mode spacing is usually the smallest energy scale, well below the relevant RG scale.The system remains integrable for a finite length ℓ, and one could use finite length form factor techniques [91,92] to investigate the effect of the finite ℓ on modes at ω ≳ ∆.
As mentioned in Section II, the UV cutoff Λ is set by the smaller of the inverse RC time of the impurity and the transmission line, and the plasma frequency of the line.A finite Λ does break integrability; however, usually, Λ ≫ T B , such that many modes, both below and above the RG scale, lie well below Λ, hence our analysis remains valid.The effect of a finite Λ could be explored by treating the cutoff terms in the Hamiltonian perturbatively, within the form factors formalism.This would require the calculation of higher order response functions -for instance, the total inelastic decay rate γ (ω) and the phase shift δ (ω) would be given by a 3-or 4-point response function (as opposed to the 2-point function considered in Section IV), which could be calculated using the formalism developed in this work.
Finally, note that while the discrete system introduced in Eq. ( 1) is not integrable due to the finite lattice spacing a, the associated scale v/a is much larger than any other scale (including Λ), and its effect on the results should be negligible.

Appendix C: S-matrix of the bulk sine-Gordon model
The S-matrix in an integrable quantum field theory is the key ingredient in the Zamolodchikov-Faddeev algebra, and reflects the purely elastic nature of the scattering in the theory.It satisfies several properties; the first is the Yang-Baxter equation, which ensures the equivalence of the factorization of nbody scattering to a product of 2-body scatterings.Unitarity and crossing symmetry imply Recall that a bar denotes charge conjugation, ± = ∓ and m = m.The S-matrix of the bulk sine-Gordon model is well-known.First, for the exchange of two (anti)solitons, where ξ = π 1/z−1 , and (C7) For integer p = 1/z, this sector of the S-matrix becomes diagonal, since S −+ +− (θ) = S +− −+ (θ) = 0, and also p S 0 (θ).Note that S 0 (θ) = −1 for the free-fermionic case z = 1/2.Next, the exchange matrix for a breather and a soliton or an antisoliton is Finally, the S-matrix for two breathers is All other terms of the S-matrix are zero.It is also useful to obtain expressions for the exchange of a right mover with a left mover in the massless limit, Sϵ ′ , as such limits of the S-matrix appear in the complete set of states in Eq. ( 21).The S-matrix is diagonal in this limit, Sϵ ′ Also note that, in the massless limit, the S-matrix of two left movers is the complex conjugate of the S-matrix of two right movers, Appendix D: Boundary reflection matrices in the Kondo and boundary sine-Gordon models The presence of the boundary introduces another component to an integrable field theory -the boundary reflection matrix, R ϵ ′ ϵ (θ), which relates an incoming state |θ⟩ ϵ with an outgoing state |−θ⟩ ϵ ′ : |θ⟩ ϵ = R ϵ ′ ϵ (θ) |−θ⟩ ϵ .First studied for solitons in Ref. [37] and later for bound states in Ref. [105], the reflection matrix can be derived from a set of axioms and properties, similar to those of the S-matrix, such as boundary unitarity: The massless limit of the reflection matrices of the bsG model was later derived in Ref. [58], where it was shown that R ϵ ′ ϵ (λ) depends only on the difference λ−λ B , where T B = e λ B is the energy scale associated with the boundary, proportional to the RG scale . The soliton reflection matrices in the bsG model are These matrices simplify considerably at z = 1/2; since R s (λ) = e λ/2 / e λ + i , we find For breathers (z < 1/2), we have The reflection matrices of the Kondo model are simpler: and It is important to note that, in the limit λ ≫ λ B , the reflection matrices become trivial in both models: for a set of excitations λ, Recall that R ϵ ′ λ ϵ λ , defined in Eq. ( 27), involves the product of scattering matrices of right and left movers given in Eq. (C10), which ensures that the phase of the product of reflection matrices in that limit is zero.In the low energy limit, λ ≪ λ B , Eq. (D8) still holds for the Kondo model, while in the bsG model we have ϵ k , (D9) where the sign factor C λ is defined in Eq. (E9).Also note that, in both models, Appendix E: Form factors in the massless sine-Gordon model

General properties
The form factors in an integrable quantum field theory can be derived from a set of axioms and conditions determined by the local conservation laws of the model, as well as the additional symmetries of the specific model.First, two excitations can be exchanged via the S-matrix:

.) . (E1)
We also have the following periodicity property: ) where s O is the spin of the operator O. Specifically, the spin of the current operators is s R = s L = 1.
The form factors are analytical functions of the rapidities in the strip 0 ≤ Imθ ≤ π, where the only singularities are two kinds of simple poles.The first are the annihilation poles, at θ i = θ j + iπ.The residue of the pole at Here C ϵ1ϵ2 is the charge conjugation matrix: where C is the charge conjugation operator.In the sine-Gordon model, m , and zero otherwise.The other residues at θ i = θ j + iπ can be found using Eq.(E1).Poles of the second kind indicate the bound states in the theory.For example, consider a form factor f O ϵ1...ϵn−1ϵn (θ 1 , . . ., θ n−1 , θ n ) in the sine-Gordon model, with ϵ n−1 = + and ϵ n = −; this form factor has a pole at θ n = θ n−1 +iθ (m) , with θ (m) = π−ξm, corresponding to a breather of type m.Its residue is proportional to the form factor f O ϵ1...ϵn−2m (θ 1 , . . ., θ n−1 , θ n ); this is known as the bootstrap axiom, which allows one to obtain form factors of breathers from those of solitons, or higher-order breathers from lower-order ones.
The crossing relations are needed to evaluate matrix elements of the form ⟨ϑ 1 , . . ., ϑ l |O|θ n , . . ., θ 1 ⟩ , which appear in the calculation of multipoint correlation functions.Following the notations of Ref. [57], we write . We also define The rapidities in the bra state of ϵ θa are analytically continued so there is no overlap between the bra and the ket rapidities, and where C θ denotes the product of the elements of the charge conjugation matrix for the excitations in the set θ: The infinitesimal imaginary part δ → 0 + removes the singularities from The equivalence of Eqs.(E7) and (E10) is guaranteed by the axioms (E1) and (E4).It is also possible to use a mixed version of the two forms, where some rapidities are analytically continued with +iδ and the others with −iδ.Choosing a specific partition ϑ = ϑ A ∪ ϑ B , we may write with the graphical representation where the sum runs over all partitions and permutations of ϑ A,B , θ.
The crossing relations are particularly simple for matrix elements of the form ⟨ ⃗ θ|O|0⟩ Specifying to the hermitian current operators, A = R, L (where s A = 1), and inserting the identity operator 1 = C −1 C between each pair of creation operators, we find where we used Eqs.(E5), (E9), C 2 θ = 1, and CAC −1 = −A.We thus have (E16)

Form factors of the current operators in the massless sine-Gordon model
The properties above hold for any integrable quantum field theory.We now focus on form factors of the current operators in the massless sine-Gordon model.These are defined by R = J 0 + J 1 , L = −J 0 + J 1 , with J µ = −ϵ µν ∂ ν ϕ.The form factors of J µ were derived by Smirnov [57], and have the general structure where g ϵ θ ⃗ θ is a function which depends only on the differences θ j − θ k .The form factors of R, L are then with ς R = −ς L = 1.The function g ϵ θ ⃗ θ vanishes when one of the rapidity differences approaches infinity, with A → ∞ and M + e A /2 → 1 in the massless limit, ⃗ θ are non-zero only if all excitations are right or left movers, respectively.This justifies keeping only the first and last rows of Eq. (21), . The form factors of R and L are related by complex conjugation and charge conjugation: If the imaginary part of a rapidity is not zero, then one has to take the complex conjugate of the rapidity as well:

.) . (E20)
A complete list of the form factors in the massive sine-Gordon model may be found in Ref. [57].Here we summarize the form factors of the right current operator R that are used in this work, taking the massless limit of the sine-Gordon model.First, the form factor of two solitons is given by where λ) , and Note that the expression for e I(λ) , which is independent of the parameter N , is valid as long as the integral converges, which depends on the imaginary part of λ (the integral expression for e I(λ) is derived from an infinite product of Gamma functions [106]).The rate of convergence may be improved by increasing N .We also have has poles in the regime 0 ≤ Im {λ 2 − λ 1 } ≤ π, corresponding to the breathers in the theory; using the bootstrap principle, we obtain the form factors for single breathers, where m is odd (f R m = 0 for even m), and θ (m) = π − ξm.If p = 1/z is an integer, S 0 has a pole at θ = iθ (m) , and then the argument of the square root has to be evaluated in the limit ξ → π/ (p − 1), In fact, one may obtain f R 1 , and also f R 111 , f R 11111 , . . .from the correspondence between the sole excitation of the sinh-Gordon model and the m = 1 breather of the sine-Gordon model.Here we only need f R 111 , The form factors above hold for any value of the coupling constant z.We also need f R +−+− , f R +−m , which are significantly more complicated for numerical evaluation.Luckily, their expressions simplify for integer p = 1/z, where The form factors become particularly simple at the free- and all other form factors are zero.
ϵ λ are at integer p = 1/z, where r 0 +−1 is evaluated.Consider Eq. ( 69).Our goal is to write the product of matrix elements, using the crossing relations in Eq. (E11), in a way that will mark the excitations according to the operators they are connected to.First, we have Here α 11 is the group of rapidities which is not passed on to the other operators, and connects the first operator does not appear in the expression above, since we simultaneously order (recall Eq. (C11)).We now multiply the above by , and use (which is the scattering matrix for left movers) and δ using Eq.(E11) and the partition ϑ A = λ 2a , ϑ B = λ 1b : Plugging everything back into Eq.( 69) using the delta functions and scattering matrices to reorder the other matrix elements, we arrive at Eq. ( 70).

Explicit expressions for the diagrams (72)-(75)
The diagrams ( 72)-( 75) correspond to the following terms: ie i(ν21+ν22−ν11−ν12)xin e i(ν21+ν22+νγ +iη)t ′′ cos (ωt ′′ ) (F7) These terms follow from the crossing relations, as shown explicitly for III in the Subsection above.It is crucial to keep track of the infinitesimals ±iδ, which determine the position of the annihilation poles of the form factors with respect to the integration contour.

The non-causal response function G cqq
ω ′ ;RRLL The non-causal response function G cqq ω ′ ;RRLL gives the response of a right-mode occupation n RR q ′ to an injected photon moving away from the boundary.It is a crucial sanity check to verify that it vanishes.The analogues of Eqs. ( 76), (77) are, in this case, Shifting κ 1 → κ 1 −iπ+3iδ in I allows to close the contour.However, now the pole at κ 1 = κ 2 − 2iη is not enclosed by the contour; hence, there are no singular contributions to the integral.This holds for II, IV as well, which means that G cqq ω ′ ;RRLL does not contribute to the spectrum, as expected.Note that the correlators comprising the response function do not necessarily vanish on their own, and we must consider their combined contribution to show that the response is zero.

Identifying and evaluating the leading contributions
Eq. ( 90) provides a general expression for the inelastic spectrum by means of a form factor expansion.The presence of breathers in the attractive regime, z < 1/2, leads to many possible terms, whose contributions are expected to decay rapidly with the number of excitations.The diagrammatic representation serves as a convenient tool to identify the leading terms.It is important to recognize that many of the form factors vanish due to the U (1) symmetry of the bulk sine-Gordon model; namely, for a mixed matrix element of the form ⟨ ⃗ ϑ|A| ⃗ θ⟩ ϵ ϑ a a ϵ θ , the total topological charge must vanish, for all m ≥ 1, and f R m1,m2 = 0 for even m 1 + m 2 , allowing us to exclude many of the terms in the expansion.The delta functions in Eq. ( 90) are also useful in the exclusion of several terms -for example, a term with ϵ 11 = ϵ γ = {} and ϵ 12 = {m} is forbidden, since the delta function δ (ω − ν 11 − ν 12 − ν γ ) implies ν 12 = ω, hence the argument of the delta function δ While one looks for terms with as little number of excitations as possible, one must be careful when mixed matrix elements of the form ⟨ ⃗ ϑ + iδ|A| ⃗ θ⟩ ϵ ϑ a a ϵ θ are involved, due to the presence of the annihilation poles; if ϑ i gets close to θ j for some i, j, the effective order of the form factor is reduced by 2 (see Eq. (E4)).That means that, in the evaluation of mixed matrix elements using Eq.(E11), it is not enough to consider only the terms which eliminate the maximal number of excitations on both sides of the matrix element (i.e.terms in Eq. (E11) with the largest number of delta functions), since the other terms in the sum are just as important.To be concrete, consider the following example: where π ± = π ± δ, and S 1+ 1+ is the S-matrix for a breather and a soliton, given in Eq. (C8).The contribution of both terms in each of the rows to some correlation function is of the same order, since and the residue of f R −+1 is proportional to f R 1 (see Eq. (E4)).The delta function δ (θ 1 − ϑ 1 ) in Eq. (F11) can be eliminated by taking the average of the two equivalent forms in Eq. (F10): Now we can expect the principal value term to be subleading with respect to the first term.
The same averaging should be applied for the diagrams drawn in Subsection V D. To illustrate this, consider the leading contribution to the z = 1/3 spectrum, depicted in Eq. (94). the above, we understand that we need to consider its "prior" diagram, from which it originates, The equivalence of the two forms is again a result of the consistency of Eq. (E11) for any choice of ϑ A , ϑ B .Now, there is a delta function term in the second diagram of each of the two forms, "hidden" within the form factor associated with b L q ′ , whose contribution cannot be neglected.The solution is the same as in Eq. (F12)average the two forms to get rid of the delta function terms.The two forms are related by flipping the order of excitations in each bra and ket states, which corresponds to taking the complex conjugate of the form factors (but not of the reflection matrices).Hence, we may get rid of the delta function by taking the real part of F in Eq. ( 82): leading to Eq. ( 95), and whose contribution to the spectrum is ) While the integration contour passes close to (but not through) the poles of f R +−+− , taking the real part of the product of form factors means that this term may be evaluated numerically as a principle value integral.Numerical inspection shows that this term is negligible compared to γ (1) 2 in Eq. ( 95).It is not always possible to get rid of all of the hidden delta function terms.For example, consider the "prior" diagrams of Eq. ( 92): and its equivalent forms can be drawn as before.We then find ) , (F18) leading to Eq. ( 93).The second and third diagrams in Eq. (F17) are mirror images of each other, and their combined contribution is where Again, taking the real part eliminates the delta function terms from f R +−+− , and the integral can be evaluated as a principle value integral; this term turns out to be negligble compared to γ (1) 1 in Eq. ( 93).However, the delta function terms cannot be eliminated from the fourth term in (F17) by taking its real part, since there are annihilation poles in both form factors related to b L † q ′ and b L q ′ .While this term involves 6 excitations and thus requires the calculation of a triple integral, it can be well approximated by considering only the delta functions, which reduce the number of excitations by 2 and therefore reduce the triple integral to a single integral:

Terms used to evaluate the spectrum
We list here all of the diagrams used in Eq. ( 90) to evaluate the spectrum.As detailed above, we take the real part of each of the diagrams, and discard all diagrams which correspond to principle value integrals, since those are subleading.The contributions indeed decay rapidly, as illustrated in Fig. 5 for the z = 1/3 spectrum.(F22) The two diagrams correspond to Eqs. ( 93) and (F20), respectively.The first diagram is the only nonvanishing term for z = 1/2, and is the only diagram used to evaluate the spectrum for z ≥ 1/2 and z = 0.47 in Fig. 3.
Appendix G: Refermionization at z = 1 /2 The bosonic Hamiltonians ( 4) and ( 8) may be mapped into Hamiltonians of free fermions at the special point z = 1/2 (known as the Toulouse point in the context of the Kondo model [107]).Solving both by means of refermionization is an important consistency check for the results of Sections IV and V.We introduce a fermionic field ψ: 2πa 0 e iπf † f e iϕ(x) , (G1) where a 0 is a short distance cutoff scale, and f † is a fermionic creation operator that anticommutes with ψ, necessary to ensure proper anticommutation relations of ψ.It is convenient to unfold the lead of length ℓ to a lead of length 2ℓ, with the impurity placed at x = 0.This allows us to expand ψ to its eigenmodes, ψ (x) = 1 √ 2ℓ k ψ k e ikx , as written in Eq. ( 50).The Hamiltonians (4) and ( 8) become quadratic under this mapping.The Kondo Hamiltonian is simpler; identifying the pseudo-spin operator in Eq. ( 8) with f , S − = f , we find where α = J 2 πa0 ℓ .The bsG Hamiltonian is slightly more complicated, as straightforward refermionization of the cosine term in Eq. ( 4) would lead to terms that are linear in ψ (x = 0) and ψ † (x = 0).To overcome this problem, one may introduce a spin operator in front of the cosine term, E J cos (ϕ (x = 0)) → S x E J cos (ϕ (x = 0)), with S x = f † + f , so that the new Hamiltonian commutes with S x and is therefore equivalent to the original Hamiltonian for either S x = 1 or S x = −1 (in the latter case, up to a shift of ϕ) [80].We then find where now α = E J 2 πa0 ℓ (we relate to the prefactors in both models as α for brevity).The mapping (G1) allows us to calculate exact bosonic correlations functions using the fermionic propagators.Working in the Keldysh formalism, we define the matrices D ab (k, k ′ ; ω), which are the temporal Fourier transform of (G7) where T is the temperature.

The total inelastic decay rate
As before, we extract the reflection coefficient from the conductance.Working in the chiral version, there is only a right moving current, and the conductance has to be calculated between x ′ < 0 and x > 0. In the Keldysh formalism, the retarded correlator is given by G (x > 0, x ′ < 0; t) = − 1 8πω ⟨T K R c (x, t) R q (x ′ , 0)⟩ .
(G8) Note the difference of a minus sign compared to Eq. ( 14), due to the unfolding of the half-infinite line.The current operator may be written in terms of the fermionic modes ψ k by inverting Eq. ( 60 (G11) This result generalizes the conductance (and hence the inelastic and elastic scattering rates) at z = 1/2 for finite temperature.At T = 0, the tanh factors become step functions, leading to This result recovers Eq. ( 51), as we identify Λ = T B /2 and Λ = 2T B for the bsG and Kondo models, respectively.The inelastic rate is again given by γ (ω) = 1 − |r (ω)| 2 , and its origin can be understood as before, this time using the language of bosonization: at z = 1/2, a photon is comprised of a fermionic particle-hole pair, both members of which scatter elastically off the boundary, but pick up different phases while doing so, leading to the splitting of the incoming photon.

The energy-resolved inelastic decay spectrum
Consider the 3-point fully retarded correlator, G cqq ω ′ (t − t ′ , t − t ′′ ) = − T K n c q ′ (t) ρ q (x in , t ′ ) ρ q (x in , t ′′ ) , (G13) where x in < 0. Note that in the chiral representation of the fermionic field, the density is defined to be rightmoving: ρ q (x in ) = 1 ℓ a=c,q k1,k2 e i(k2−k1)xin ψ a † k1 ψ ā k2 . (G14) The bosonic occupation number n q ′ is written in terms of the fermionic operators as n c q ′ =b q † q ′ b q q ′ + b c † q ′ b c q ′ = b q † q ′ b q q ′ + π ω ′ ℓ k1,k2 a1,a2=c,q ψ a1 † k1+q ′ ψ a1 k1 ψ a2 † k2 ψ a2 k2+q ′ .(G15) Note that, plugging the above to Eq. (G13), the b q † q ′ b q q ′ term leads to an all-quantum Keldysh correlator, which vanishes identically.We thus have   where we sum over all permutations P of {1, 2, 3, 4} (Pi denotes the permutation value for the index i), and ζ P = 1 (−1) for an even (odd) permutation.Here we only consider connected terms, where each ρ leg is connected to both b † q ′ and b q ′ legs -namely, |P1 − P2| > 1 or |P3 − P4| > 1 (one may show that the disconnected terms are proportional to δ (ω − ω ′ ) and therefore correspond to elastic scattering, similarly to Eq. ( 91)).Note the subtraction of the background term in Eq. (G22), similarly to the subtraction of the Kronecker deltas in Eq. (83).Evaluating the sums over k 1,2 , taking care to close the integration contours in the half planes allowed by the x in exponentials, and keeping only contributions that are singular in η, we find We can now plug G cqq ω ′ into Eq.( 54), thus obtaining the spectrum at z = 1/2 at finite temperatures.At T = 0, the tanh functions become step functions; direct inspection of the sums over the permutations shows that which are both identical to the form factors results, and are simple enough to lead to closed analytical expressions.Indeed, one may consider Eq. ( 93), and use the form factor in Eq. (E32) and the reflection matrices in Eqs.(D3) or (D6) to recover the expressions above.

Figure 3 .
Figure 3. (a) The energy-resolved inelastic decay spectrum as function of ω ′ , at several fixed values of ω/TB, for the bsG and Kondo models and z = 1/3, 1/2.The diagrams used to evaluate the spectrum are listed in Appendix F. (b) The ratio between the LHS and RHS of Eq. (53) for both models and several values of z.Note that the power laws of γ (ω) are recovered by the sum rule for all z.
This form allows us to decouple the inductive coupling between the array and the fluxonium by applying a unitary transformation, U f = e iα φf Qtot , with Qtot = N n=0 Qn , which shifts the array phases, φn → φn + α φf for all 0 ≤ n ≤ N , as well as the fluxonium charge, Qf → Qf − α Qtot .The array-fluxonium coupling becomes capacitive, and the Hamiltonian reads H = Hf + Ha + Hc , with Hf

3 e
) Again, F (λ) does not depend on N , which is a useful parameter to increase the rate of convergence.The pole of f R 111 at λ 3 − λ 2 = iξ yields the form factor f R 12 , λ1 + 2 cos ξ 2 e λ2 e λ1 e λ2 e 2λ1 + e 2λ2 + 2 cos ξ 2 e λ1 e λ2 × F λ

Figure 5 .
Figure 5.The weight of the contributions to the z = 1/3 spectrum, evaluated using the sum rule.