Three-body bound states in antiferromagnetic spin ladders

Stable bound quantum states are ubiquitous in nature. Mostly, they result from the interaction of only pairs of particles, so called two-body interactions, even when large complex many-particle structures are formed. We show that three-particle bound states occur in a generic, experimentally accessible solid state system: antiferromagnetic spin ladders, related to high-temperature superconductors. This binding is induced by genuine three-particle interactions; without them there is no bound state. We compute the dynamic exchange structure factor required for the experimental detection of the predicted state by resonant inelastic X-ray scattering for realistic material parameters. Our work enables us to quantify these elusive interactions and unambiguously establishes their effect on the dynamics of the quantum many-particle state.


Introduction
Binding phenomena were at the origin of the development of quantum mechanics. The quest to understand why electrons flying around positive nuclei form stable entities, namely atoms, has led to a completely new foundation of physics over hundred years ago. The formation of the structures from nuclei and electrons over atoms, small and large molecules to macroscopic liquids and solids is founded on quantum mechanical binding. But also down to smaller and smaller length scales binding is crucial in the formation of the atomic nuclei formed by protons and neutrons due to the strong force and of heavy particles in the Standard Model of high-energy physics leading to hadrons, namely mesons, bound states of two elementary quarks, and baryons, bound states of three elementary quarks.
Binding effects are especially dominant in correlated antiferromagnetic spin system in low dimensions. Such systems have already revealed some fascinating binding phenomena before. A chain of localized spins S = 1/2 coupled on adjacent sites and subjected to a magnetic field displays fractional excitations ((anti-)psinons) of halfinteger quantum numbers, spin flips (magnons) with integer spin number, and in particular bound states of n magnons (n-strings). These states were theoretically predicted up to 90 years ago [1][2][3][4] , but have been verified experimentally only recently 5,6 . Spin ladders, in particular, are solid state systems which are intensively studied due to their relation to high-temperature superconductors 7 . Moreover, they represent a paradigm for an entangled magnetic many-body systems without long-range order 8 . Two-body bound states in spin ladders have been measured for S = 0 mediated by a phonon in infrared absorption 9 and for S = 1 by inelastic neutron scattering 10 .
In spite of the abundance and complexity of bound states most of the binding phenomena are induced by two-body interactions. This means that the underlying potential only depends on the positions and properties of pairs of particles. This is particularly obvious in chemistry where all atoms and molecules are held together essentially by Coulomb potentials. It is the paradigm of a two-body interaction which can bind not only two, but almost infinitely many particles.
In this communication, we show that there are also existing, realistic solid state systems where a genuine (irreducible) three-body interaction provides the vital extra for the formation of a bound state. This bound state would not exist if only two-body interactions were present. The key point is that an irreducible three-body interaction acts only if three particles are present. It has no effect on two particles or a single particle. We explicitly calculate the effect of this bound state for dynamical correlations in spin ladders using continuous unitary transformations (CUT). We find an effect stronger than 50% for experimentally relevant observables that can be measured with resonant inelastic X-ray scattering (RIXS).

Theoretical model
We consider the antiferromagnetic spin ladder 7,8 where S i,τ denotes the vector spin operator at rung i along the leg and on site τ = 0, 1, see Fig. 1(a). Spin ladders with x := J leg /J rung 1 are realized in cuprates to a high degree of accuracy 7 so that experimental verification of our predictions is possible. For x = 0, the excitations are local S = 1 triplets above the S = 0 singlet ground state ( Fig. 1(b)). For x > 0, the elementary excitation is no longer localized on one rung only, but it is smeared out over a number of rungs, the size of which is given by the correlation length ( Fig. 1(c)). It is now called a triplon [11][12][13][14][15] . Note that triplon excitations are typical also in other dimerized systems, such as in Shastry-Sutherland magnets 16 . In contrast to the number of triplets, the number of triplons is conserved also for finite x and therefore, triplons provide a natural basis for the description of dimerized systems 17,18 . For clarity, we distinguish in this article between the initial non-conserved triplets and the conserved triplons. In literature, this distinction is not always used [17][18][19] .
We stress that the original triplets are not the appropriate elementary excitations because already in the ground state there is an infinite number of them admixed. Hence, it does not make sense to refer to a one-, two-or three-triplet state since any eigen state comprises an infinite number of them. Instead, we have to use the elementary excitation resulting from renormalization. This is analogous to the vacuum fluctuations in quantum field theories. The observed and measurable elementary particles are quasiparticles, dressed by vacuum fluctuations, whose properties stem from renormalization.
The leg coupling also leads to an attractive interaction between pairs of triplons which are not far apart from each other, see Fig. 1(d). The competition between attractive interaction and kinetic energy determines the formation of bound states in the spin ladders. As triplons are equivalent to mobile S = 1 spins that interact antiferromagnetically on the spin ladder, it is expected that the energetically lowest bound state is in the total S = 0 sector. Two-body bound states with total spin S = 0 and S = 1 induced by the two-triplon attraction are established theoretically 11,12,14 and have been measured in experiments 9,10 . The aim of the present communication is to show that three-triplon bound states occur, see Fig. 1(e), which are essentially due to genuine three-triplon interactions.
The occurrence of the three-triplon interactions results from a mechanism very analogous to the attractive interaction between electrons resulting from the exchange of phonons which eventually leads to the formation of Cooper pairs and superconductivity. The propagation of an electron excites a phonon which is captured by another propagating electron inducing their attraction 20,21 . Similarly, Fig. 2(b) illustrates that the creation of a pair of triplons, triplon propagation and subsequent pair annihilation induces a genuine three-triplon interaction. Besides pair creation and annihilation, i.e., the vacuum fluctuations, the hardcore property of triplons, excluding more than one per dimer, generates the interaction involving three triplons on three dimers. Otherwise the interaction would be a single-particle irreducible one, see Fig. 2(a). We stress that this mechanism is not specific to one dimension, but applies in any dimension to systems of coupled dimers 16,22,23 and, even more generally, to any system with finite Hilbert space dimension at each site because the latter implies a hardcore property, i.e., for spin flips.
To obtain a quantitatively correct description of the elementary excitations in the spin ladder we express the original model (1) in terms of triplet creation and annihilation operators, see Methods. Then, by a systematically controlled change of basis, the model can be mapped to one in which the number of triplons is conserved. This is a vital step to make further analysis feasible and it is achieved by means of a CUT 21,[24][25][26][27] . The CUT is defined by the infinitesimal generator η which is classified by a label indicating whether the 2-triplon subspace is separated from all subspaces with n > 2 triplons (label (2:n)) or even the 3-triplon subspace (label (3:n)) is separated as well. The resulting set of differential equations describes the renormalization flow of the dispersion and the interactions. These differential equations are truncated in some order in x to keep their number finite and finally solved numerically. In the course of the CUT the interactions of the triplons are renormalized and novel types of interactions are generated. In particular, the two-triplon interaction is modified and three-triplon interactions are induced. This is the mechanism how a genuine, irreducible three-body interaction comes into play. The possibility to distinguish the influence of the genuine two-triplon interactions from the one of the genuine three-triplon interactions is the crucial ace of the CUT approach. An extensive derivation of the applied methods can be found in the Supplementary Methods.

Theoretical description of RIXS response
To identify the effect of three-triplon interactions, the symmetries of spin ladders help greatly: reflection about the centerline defines even and odd parity and the spin isotropic model conserves the total spin. Since triplons are of odd parity 11,12,14 the excitations in the odd channel consist of one or three or five triplons and so on. If in addition one uses a spin conserving (SC) probe the total spin of the excited states is zero as in the ground state. A triplon has spin S = 1, hence in the odd S = 0 channel the leading contribution is the one formed by three triplons only. This channel can be addressed by light scattering, for example RIXS 28 , or absorption in the terahertz (THz) range 5,6 . RIXS has the advantage of substantial momentum transfer due to the high energy of the scattered photons, in particular for hard X-rays as used for the Cu K-edge, while THz absorption provides high energy resolution, but stays at zero momentum. RIXS spectra can be resolved into spin-conserving (∆S = 0) and spin non-conserving (∆S = 0) at the Cu L 3 -edge of cuprates 29,30 . The spin-conserving channel can be well captured by the dynamical exchange structure factor (DESF). In contrast to the Cu L-edge, the channels at the Cu K-edge and oxygen K-edge in cuprates are purely spin conserving 31,32 , which is advantageous for experimental verifications of our predictions. RIXS resolution for hard X-rays has also improved considerably in the recent past 33,34 .
To compute the RIXS response at the Cu Kand L-edge of cuprates, we use the dynamical correlation functions given by the ultra-fast core-hole liftime (UCL) approximation 35,36 . We focus on the SC channel, not accessible by inelastic neutron scattering which measures the dynamics structure factor (DSF). For data of the DSF see Supplementary Note 1. The SC channel is captured by the DESF given by is the spin exchange observable, |g and | f are the ground and final states with energies E g and E f , respectively, and ω is the energy loss to the system. Usually, the DESF is evaluated by exact diagonalization (ED) 35,[37][38][39] or by density matrix renormalization group 40 . These approaches provide spectra consisting of a multitude of finite-size peaks. But they do not allow to trace back the physical origin of the spectral features. The distinction between continuous scattering contributions and peaks from bound states can also be challenging.
Scattering states of two or more triplons lead to a continuous contribution to the DESF at fixed total momentum transfer q. This holds also for a two-triplon bound state scattering with a single triplon. Only a bound three-triplon state yields an infinitely sharp δ -peak in S ex (q, ω) at given q. Thus this is the smoking gun feature we have to look for. For rendering purposes we will broaden it slightly; but in experiments it will show up as sharp peak limited only by the resolution of the apparatus.

Three-triplon bound states
We provide data for x = 1.2 and x = 2 because these values represent the experimentally relevant range in cuprates 9, 38, 41 and the telephone number ladder La 5.2 Ca 8.8 Cu 24 O 41 in particular. Figure 3 depicts the DESF at (q x , q y ) = (π/(2a), π/a). Panels (a) and (b) show the three-triplon bound states that appear separated below the three-triplon continuum for x = 1.2 and x = 2.0. The DESF of the bound state is plotted in red and the continuous 3/15 DESF in orange. The latter stems from three asymptotically free triplons and from scattering of a S = 1 two-triplon bound state with a single triplon.
We highlight that the three-triplon bound states have significant weight compared to the continua: for x = 2.0 over 50% of the spectral weight resides in the three-triplon bound state. The energy separation between the bound state and continuum is small for x = 1.2 while it is becomes substantial for x = 2. This makes compounds with higher ratios x interesting for the experimental verification of bound three-triplon states. For reference, CaCu 2 O 3 is known to be a host for a spin ladder with large ratio x 42 .
To explicitly show the effect of the genuine three-triplon interactions panels (c) and (d) of Fig. 3 display the DESF but without the three-triplon interactions. Then, only a square-root divergence of the DESF appears at the lower band-edge. For x = 2.0 there may exist an extremely weakly bound state. But the distinctive three-triplon bound states, clearly separated from the continuum only occurs as a consequence of the irreducible three-triplon interactions. This is the key result of our investigation. Figure 4 depicts the DESF at q y = π/a as function of q x . Because of the symmetry S(q x , q y , ω) = S(−q x , q y , ω) it is sufficient to show the positive half of the Brillouin zone. In all cases, the total weight is largest around q x = π/(2a). While Fig. 4 (a) and (b) depict the full calculation with all interactions, Fig. 4 (c) and (d) exclude the irreducible three-triplon interactions. The solid white line shows the lower edge of the three-triplon scattering continuum if only the one-triplon dispersion is considered, i.e., three asymptotically free triplons are involved. The dashed white lines depict the actual lower continuum edge including scattering states from an S = 1 two-triplon bound state and a single triplon. Hence the difference between the solid and the dashed line for q x < 0.65π (x = 1.2) and q x < 0.75π (x = 2.0), respectively, stems from the irreducible two-triplon interactions. We emphasize that only binding can induce states at lower energies. Without binding, interactions can shift spectral weight, but only between the edges of the continua. This is indeed one effect of the irreducible three-triplon interactions: they shift the weight significantly towards the lower band edge as can be clearly discerned by comparing Fig. 4 (a) with Fig. 4 (c) and Fig. 4 (b) with Fig. 4 (d).
The energy of the bound state formed from three triplons lies below the lower continuum edge for q x < 0.65π/a (x = 1.2) and q x < 0.75π/a (x = 2.0), respectively. These bound states are close to the lower continuum edges as expected from Fig. 3. Their dependences on q x have a similar shape. Figure 3 indicated that the spectral weight of the bound state is significant relative to the weight in the adjacent continuum. This feature holds generally for most values of q x as shown comprehensively in Fig. 5 for four different values of x as function of q x . The ratio of the spectral weight of the bound state to the spectral weight of all scattering states I bound /I cont is plotted. We stress that the weight of the scattering states is integrated up to the highest energies far beyond the energy of the three-triplon bound state. The maximum of the relative spectral weight moves to higher q x upon increasing x, i.e., increasing J leg . Also the relative weight increases with increasing x. Even for low values of x, the maximum of the weight of the bound state is never an order of magnitude below the weight in the continuum. For x = 2.0, the bound weights even exceed the continuum weights.
For completeness, the Supplementary Note 1 provides further data for the DESF and Supplementary Note 2 provides an estimate of the importance of the contributions of four triplons and a comparison to data for finite spin ladders obtained by ED. Comparison to the sum rule obtained from ED also shows, that the three-triplon sector contains almost all of the spectral weight for all x considered here. Additionally, data for the standard dynamic structure factor is shown which is relevant for inelastic neutron scattering.

Discussion
Binding effects are fundamental to understand the structure of the surrounding matter ranging in length scale from femtometers to meters. Most bound states stem from two-body interactions such as the Coulomb potential. But in the present study, we have shown that in a realistic, correlated condensed matter system genuine (irreducible) three-body interactions between the elementary excitations are crucial as well. This is achieved for a generic antiferromagnetic spin ladder as it is realized in cuprates well-known from the field of high-temperature superconductivity, but the underlying mechanism applies to wide classes of lattice systems in any dimension.
We systematically derived the three-body interactions of the elementary triplons by CUTs. We identified a probe channel (odd parity and S = 0) in which the three-triplon states represent the leading contribution. We established that this channel can be probed by RIXS and computed the DESF relevant for RIXS. A particular asset of the CUT approach is that one can switch on or off the genuine three-body interactions. In this way, we showed that it is the three-triplon interactions which induce a significant shift of spectral weight in the DESF to lower energies. Most notably, a bound state formed from three triplons appears in a large part of the Brillouin zone -but only if the three-triplon interaction is taken into account. Although it is close in energy to three-triplon scattering states its weight is significant and partly dominates over the weight of the scattering states.
We highlight the fundamental difference of our finding to completely frustrating, antiferromagnetic diagonal couplings in spin ladders, where the number of triplets is already conserved so that triplets and triplons coincide. The complete frustration leads to an enhanced attraction of adjacent pairs of triplons such that three-triplon bound states (Fig. 1(e)) and even states of n-triplon strings (Fig. 1(f)) are predicted 43 . These local bound states stem from two-triplon interactions; no three-body interactions are involved. Our key observation for unfrustrated spin ladders is that even for small x the three-triplon interaction is sufficiently strong for the formation of a three-triplon bound state with small energy separation from the continuum. For x 1 this state gains significant spectral weight with a clear energy separation.
These results suggest that an experimental detection of the effects of three-body interactions is possible. Nextgeneration RIXS facilities with improved energy resolution should be able to tackle this challenge.  10 , with the strongest response at q x = π/(2a) in the Brillouin zone. Alternatively, we point out that verification is also conceivable by THz absorption which has a very high resolution in energy. Here the challenge is that THz light hardly has any momentumhk so that only states with zero momentum are probed where no significant effects of the three-triplon interactions appear. But if a spin ladder material showed a slight distortion with periodicity of four rungs, the states at q x = mπ/(4a) with m ∈ {1, 2, 3} are folded to the center of the Brillouin zone and three-triplon binding would become detectable 44,45 , as for instance shown for screwed spin chains 5 . In these systems, the oxygen K-edge, which probes spin-conserving channel exclusively, can also be used in spite of it having access to restricted momentum transfer 32,39 .
The theoretical results open up the stage for the study of the effects of genuine three-body interactions occurring in generic lattice models for condensed matter. In one dimension in particular, where quantum effects are typically strongest, they could be an important ingredient for the formation of n-triplon strings ( Fig. 1(f)) generalizing the previously detected Bethe strings 5, 6 to a much broader class of solid state systems not restricted to integrable systems such as spin chains. Our results open perspectives for ways to investigate complex bound states, connecting to various cross-boundary efforts with the same aim, such as the discovery of three-body correlations in ultracold atoms 46 , the formation of topological bound states in non-Hermitian systems of photonic lattices 47 and bound states for quantum computation in superconducting nanowires 48 . With recent improvements in experimental techniques, we are hopeful of verification of our findings in the near future.

Methods
Further details on the methods and numerical details can be found in the Supplementary Methods.

Triplet and Triplon representation
Starting from the regime of strong rung couplings, we reformulate (1) in terms of triplet creation operators t α, † i |s = |t α i where |s is the singlet ground state for x = 0 and |t α i is an S = 1 triplet state of flavor α ∈ {x, y, z} on rung i. The triplet annihilation operator t α i annihilates the specific triplet at rung i and restores the singlet on this rung. The triplet operators satisfy the hardcore boson algebra t α i ,t . Representing the spin operators by the triplet operators allows us to map the spin systems onto an interacting quasiparticle problem 21 . In the course of the CUT, the triplet operators are gradually converted to conserved triplon operators. Note that we distinguish between unconserved triplets and conserved triplons here. This distinction is not

5/15
always utilized in literature. Explicit expressions for the Hamiltonian and the observables in the triplet language are given in the Supplementary Methods.

Foundation of CUT
The basic idea of a dedicated unitary transformation is to change the basis such that the problem under study becomes easier to tackle. The continuous transformation has the advantage to perform the basis change in a renormalizing fashion, i.e., processes linking large energy difference are transformed most quickly and processes with smaller and smaller energy differences are eliminated slower and slower in the flow of the renormalization parameter → ∞. This is important because one can never track the complete Hamiltonian along the unitary flow, but truncations are necessary. The renormalizing property helps to obtain robust effective models 21,24,26,[49][50][51] .
We use CUTs here to eliminate the terms in the Hamiltonian which change the number of particles, i.e., the number of triplons. In this way, this number becomes a conserved quantity, the ground state is the vacuum of triplons and hence the calculations of dynamic zero-temperature correlations are simplified to tractable few-particle calculations. To evaluate the spin conserving contributions to the RIXS spectra we decouple the interacting multitriplon sectors from each other by applying the CUT to the Hamiltonian and to the relevant observables. The resulting effective Hamiltonian consists of n-particle irreducible parts, i.e., the one-triplon dispersion, the twotriplon interaction, and the three-triplon interactions which are of particular interest here. Irreducible higher triplon interactions, e.g., of four triplons, also occur, but only matter if four or more triplons are considered.

CUT generator scheme
The CUT provides a set of differential equations for the prefactors of a large number of interactions which proliferate exponentially with system size. Hence truncations are required which systematically control the accuracy of the obtained effective Hamiltonian and effective observables. Here we employ two concepts. First, we do not separate all m-triplon states from all n-triplon states with m = n, but use either the (2:n) or the (3:n) generator. The generator (m:n) separates the 0, 1, . . . m triplon states from all higher states with n > m 27 , i.e., we decouple either the zero, one, and two triplon states from the rest or we decouple the zero, one, two, and three triplon states from the remaining states. The latter, however, does not work up to very high order in x, see below.
Second, we use the expansion parameter x = J leg /J rung to classify the various processes. We keep all processes which influence the effective Hamiltonian up to a certain order in x. The physics behind this idea is that a larger order translates to a larger range in space, i.e., to a larger distance between the interacting triplons on the rungs. This concept leads to large, but tractable sets of differential equations which then are integrated numerically. The whole procedure is dubbed directly evaluated enhanced perturbative CUT or deepCUT for short and has proven to be reliable and efficient 21 up to x = 3. The (3:n)-generator was used to calculate the irreducible three-triplon contributions. If the highest robust order was smaller than 10, the (2:n)-generator was used to calculate the irreducible two-triplon and one-triplon contributions in order 10.
The accuracy of the CUT approach is examined by comparing it to lower order CUT calculations as well as to results from ED on finite systems in the Supplementary Note 2.

Details on dynamic structure factors
Based on the effective observables obtained from the CUT we use Lanczos tridiagonalization to compute the continued fraction expansion of the DSF and the DESF, see Supplementary Methods. This renders finite-size effects negligible. At the same time we can identify the physical processes that contribute to the spectral functions. A square-root terminator was employed to terminate the continued fraction so that the contribution of the scattering states is found as continuous function without resorting to artificial broadening 52 .
The bound states are mathematically δ -distributions in the spectral densities, i.e., they are infinitely sharp, which prevents a graphical representation. Thus, the bound states (and only them) are artifically broadened in our plots. They are displayed as Lorentzians of width 5 · 10 −4 J rung .
The lower edge of the three-triplon scattering states including two-triplon interactions is calculated from the converging values of the Lanczos coefficients 52 . The lower edge at wave number q without any interactions is  Origin of three-triplon interactions. The term is depicted in real space at dimer r and interdimer distances δ , δ , δ ; note that this term arises in any dimension and for any lattice model with finite dimensional local degrees of freedom. Finite x implies hopping, pair creation and annihilation processes during the renormalization by CUT. The blue arrows indicate the incoming triplons, red the scattered triplons and the black arrows internal triplon propagation. For normal bosons (a), the combined process is single-particle irreducible and corresponds to an effective hopping. For triplons (b), the hardcore constraint (black circles) induces three-triplon interactions in leading order x 3 . The explicit expression for δ = 1, δ = 2, δ = 3 is given in the Supplementary Methods.  Ratio of bound to continuum excitations weight. The ratio is computed using the spectral weight I bound of the three-triplon bound state and the weight I cont of the continuum of three-triplon scattering states for q y = π/a. The groups of ellipses with the letter 'T' denote three-triplon states. Asymptotically scattering triplons are colored orange and bound triplons are colored red.