Phenomenology of Relativistic $\mathbf{3} \to \mathbf{3}$ Reaction Amplitudes within the Isobar Approximation

Further progress in hadron spectroscopy necessitates the phenomenological description of three particle reactions. We consider the isobar approximation, where the connected part of the $\mathbf{3}\to\mathbf{3}$ amplitude is first expressed as a sum over initial and final pairs, and then expanded into a truncated partial wave series. The resulting unitarity equation is automatically fulfilled by the $B$-matrix solution, which is an integral equation for the partial wave amplitudes, analogous to the $K$-matrix parameterization used to describe $\mathbf{2}\to\mathbf{2}$ amplitudes. We study the one particle exchange and how its analytic structure impacts rescattering solutions such as the triangle diagram. The analytic structure is compared to other parameterizations discussed in the literature. We briefly discuss the analogies with recent formalisms for extracting $\mathbf{3}\to\mathbf{3}$ scattering amplitudes in lattice QCD.


Introduction
Modern high-energy experiments are accumulating high quality data on three-hadron final states, that are expected to be the main decay channels of several poorly known or missing resonances. These include, for example, the enigmatic a 1 , ω 2 , and the exotic π 1 resonances that can be studied in peripheral production at COMPASS, GlueX, and CLAS12 [1][2][3][4][5][6][7]. In addition to conventional hadrons, many of the exotic XY Z and pentaquark states observed in the heavy quarkonium sector [8][9][10], are also found in three particle final states. a e-mail: ajackura@iu.edu Many of these newly observed or anticipated states lie close to thresholds of their decay products. For example, the mass difference between the X (3872) [11] and the D 0D0 π 0 threshold is only 6 MeV. The proximity of the three particle threshold together with the possibility of long-range interactions mediated by a single pion exchange can significantly influence the X (3872) line-shape [12] and one needs to carefully analyze the role of pion exchange and whether it is able to bind D * 0 andD 0 [13][14][15][16][17]. In the light meson sector, the recently observed a 1 (1420) [2] is yet another candidate for a state not expected in the quark model that can be influenced by meson exchange interactions and thresholds [18,19].
The problem of constraining 3 → 3 reactions from the Smatrix principles of unitarity and analyticity has been studied previously in Refs. [40][41][42][43][44][45]. In this paper we extend these earlier works and clarify some of the results. Moreover, we present the 3 → 3 reaction amplitudes in a way that can be directly translated to the finite volume.
Our description relies on the isobar approximation, where the amplitude is constructed as a sum of truncated partial wave expansions. This provides a good description of threeparticle final states in the resonance region, where analyses of Dalitz plots indicate that they are dominated by intermediate two-body resonances. For example, the decay of the a 1 (1260) resonance into three pions occurs primarily via a decay to the ρπ intermediate state with the subsequent decay of ρ to two pions [3,4]. The isobar approximation can be regarded as an effective way to incorporate the relevant singularities in all Mandelstam variables, and will be discussed in detail later.
The rest of the paper is organized as follows. In Sect. 2 we define the 3 → 3 amplitude for three spinless particles and discuss the relevant kinematics. In Sect. 3 we introduce the isobar approximation and investigate the consequences of unitarity. We explain the difference between isobar and the partial wave amplitudes , which are often confused. In short, we use the isobar representation to describe the 3 → 3 amplitude, A = A k j , where the indices k and j label the spectator particle in the final and initial state, respectively. We refer to the A k j 's as isobar-spectator amplitudes, since they can be pictured as scattering of a quasi-particle, the isobar, and a stable spectator. The latter are expanded in partial waves of the three-particle system. Unitarity constrains the 3 → 3 amplitudes on the real energy axis, which results in specific relations involving the imaginary parts of the partialwave-projected isobar-spectator amplitudes. Unitarity alone does not uniquely specify partial wave amplitudes, as evident, for example, in the K -matrix parametrization of 2 → 2 scattering amplitudes [46,47]. In Sect. 4 we discuss a specific parameterization for the isobar-spectator amplitudes which satisfies the three-body and two-body unitarity. It is given as a solution of a set of linear integral equations that involve, among others, the one particle exchange (OPE) as a driving term. We call this the B-matrix parameterization and it satisfies, where B is the driving term that contains the OPE, τ is a known function of the phase space and of the 2 → 2 amplitudes. The product formally represents an integration over the intermediate isobar mass. In contrast to Ref. [45], we restrict the domain of the integrals to physical values of energies. This enables us to use the experimentally accessible subchannel amplitudes and we also discuss the consequences of this restriction. We derive Eq. (1) for isobars with arbitrary spin s, and for any value of the isobar-spectator orbital angular momentum . The B-matrix parameterization can be analytically continued to the complex energy plane and in Sect. 5 we discuss aspects of its analytic properties. Specifically, the one parti-cle exchange process has some unique features, as it contains a kinematic singularity due to the exchange of a real particle, which can be isolated from the full 3 → 3 scattering amplitude. In addition, we also study the triangle amplitude that emerges from the B-matrix parameterization, and the relation to the Bethe-Salpeter solution. We summarize our results in Sect. 6.

Kinematics, invariants, and amplitudes
We consider elastic scattering of three distinguishable, spinless particle, e.g. DπD, K πK , or π + π − π 0 . The particles have mass m j , where j = 1, 2, or 3 labels the individual particles. A single particle state, with four-momentum p j = (ω j , p j ), where ω j = m 2 j + |p j | 2 is the energy and p j is the three-momentum, is denoted |p j and has relativistic normalization p k |p j = (2π) 3 2ω j δ (3) (p k − p j )δ k j . We are interested in the S-matrix element of the elastic 3 → 3 scattering process. We can decompose the S-matrix as S = 1 + i T . The T -matrix contains two terms, T = T d + T c , where the disconnected part, T d , involves interactions of two particles at the time with the third one being a spectator, while the connected part, T c , contains interactions of all three particles. The disconnected part can always be identified kinematically by the spectator momentum conserving delta function [48]. The disconnected part is written as T d = j 1 j ⊗ T ( j) , where 1 j is the identity operator in the single particle space of the spectator, j and T ( j) describes 2 → 2 scattering between the other two particles. The amplitudes associated with the matrix elements of scattering operators T d and T c are defined as F and A, respectively. Specifically, the connected amplitude A is given by p | T c |p = (2π) 4 δ (4) (P − P)A(p ; p), (2) where |p ≡ |p 1 p 2 p 3 and |p ≡ |p 1 p 2 p 3 denote the initial and final states of the three particles, and P = p 1 + p 2 + p 3 and P = p 1 + p 2 + p 3 are the initial and final total four-momenta, respectively, as illustrated in Fig. 1. Time-Reversal symmetry implies that the amplitude is symmetric in the initial-final state variables, A(p ; p) = A(p; p ). The chosen normalization implies that the amplitude A(p ; p) has mass dimension −2. The disconnected amplitudes F j are defined by where the delta function enforces that the spectator j does not interact. We remark that the F j is the genuine 2 → 2 scat- tering amplitude, as required by the LSZ construction [48]. We also define P j ≡ P − p j and P j ≡ P − p j as the initial and final total four-momenta of the interacting pair recoiling against spectator j, cf. Fig. 2.
In this paper we adopt the so-called spectator notation or odd-one-out notation [49], where the 2 → 2 amplitudes associated with the spectator j are labeled by the spectator index. The spectator notation requires additional information specifying the first particle in the two-particle system. There are two conventions which are useful for our discussions: the two-pair convention, and the cyclic convention. The twopair convention is more practical when interaction in one of the three pairs is negligible. An example of such a system is π + π + π − , where the π + π + system interacts weakly. In this case it is convenient to choose the noninteracting system as, say, particles (13) and designate particle 2 as the second particle for both the interacting sub-systems. Therefore, the spectator index j = 1 and j = 3 uniquely identifies the two orderings in the pairs to be (32) and (12), respectively. If the interactions in all three subchannels are important, one can define the ordering through cyclical permutation, i.e. the spectator label j = 1, 2, 3 corresponds to ordering of the two particles subsystems as (23), (31), and (12), respectively. For simplicity, in the following we assume only two relevant subchannels, and use the former convention. Generalization to the latter case is straightforward. The type of applications we have in mind are systems like MMπ elastic scattering, where M is an open-flavor meson, such as K , D, and B. The interacting pairs will be assumed in the Mπ andMπ channels only, and pion being designated as particle j = 2.
The 3 → 3 amplitude depends on eight independent kinematic variables. The choice of variables largely depends on the kinematical range of interest, e.g. the low vs high total energy region. Here we are interested in the low-energy region and use the following redundant set of Mandelstam variables, where s, σ j , and σ k are the invariant mass squares of the total three particle system, the initial pair, and the final pair, respectively. The transferred momenta, t jk and u jk , are between the initial and final spectators and the initial pair and final spectator, respectively. The Mandelstam invariants are related by energy-momentum conservation, In the physical region of the 3 → 3 reaction, s can take any value above the three particle threshold, s ≥ s th = (m 1 + m 2 + m 3 ) 2 , while the subchannel invariant masses σ j and σ k are bounded by σ are the sub-energy thresholds, e.g. σ (th) We will need the relations between the invariants and energies and scattering angles, in three reference frames. The frames of interest will be the overall center-of-momentum frame (CMF) and the isobar rest frame (IRF). There are two IRFs corresponding to the initial and final states: the initial IRF j , labeled by the spectator j, and the final IRF k , labeled with spectator k and a prime.
To distinguish momenta in the CMF we denote them by a , i.e. P = P = 0. In the CMF the scattering angle, k j , is defined as the angle between the initial and final state spectator momenta, cos k j ≡ p k · p j , where p j and p k denote the CMF orientations of the initial and final spectators, respectively. The kinematic variables in the other frames, IRF j (P j = 0) and IRF k (P k = 0) are obtained from the CMF by a Lorentz boost in the direction opposite to momentum of the corresponding spectator. The momentum of the first particle in the pair is denoted by q j , and q k in IRF j and IRF k , respectively. Orientation of these momenta are given by solid angles, q j = (γ j , χ j ) and q k = (γ k , χ k ), respectively. Here, γ j and γ k are the azimuthal angles between the decay plane of the isobar and the isobar-spectator scattering plane and χ j and χ k are helicity angles, cf. Fig. 3 for the specific scheme (23)1 → (12)3. The relations between all Fig. 3 Connection between the three reference frames for the (32)1 → (12)3 system. The total reaction plane in the CMF is shown in gray, and the initial and final IRFs are shown in blue and green, respectively. The Lorentz boosts β 1 and β 3 indicate the transformations between the three frames relevant kinematical variables and the Mandelstam invariants are given in "Appendix A". In the following, we will use the set ( q k , σ k , s, t jk , σ j , q j ) to describe the isobar-spectator amplitude.

Unitarity relations
We consider elastic unitarity in the physical region of the 3 → 3 reaction below inelastic thresholds. It yields two relations [50], one for the disconnected 2 → 2 amplitude F j and one for the connected 3 → 3 amplitude A. For F j , one finds where is the phase space for the two particle system, and q j is the intermediate state relative momentum. The IRFs are defined with their z-axes defined along the opposite direction of the spectator and their x-axes defined by their azimuthal angles w.r.t. the total CMF plane spanned by the initial and final spectator momenta, cf. Fig. 3. Note that from energymomentum conservation, |q j | = |q j | = |q j |. Figure 4 is a diagrammatic representations of Eq. (6). Elastic unitarity yields the following condition for the connected 3 → 3 amplitude, where μ jk is the mass of the exchanged particle that is neither j nor k, e.g. if j = 1, and k = 3, then the exchanged mass is μ 13 = m 2 . Note that the evaluations p k = p k in the second and fourth lines enforce that σ k = σ k , and similarly in lines three and four, p j = p j implies that σ j = σ j . Figure 5 is a diagrammatic representation of Eq. (8) and its derivation is given in "Appendix B". The implications of unitarity for the F j are summarized below. The unitarity relation for the connected, A amplitude is more complicated. The first term in Eq. (8) is analogous to the 2 → 2 case, in the sense that it is given by the product of the same connected amplitude A. The next two terms originate from the contribution to S † S given by the product of T c and T d , and represents the situation when only two of the three particles rescatter. The last term is the contribution to the imaginary part of the connected amplitude from the product of two disconnect amplitudes and reflects the real one particle exchange process. Since the unitarity relation deals with physical, on-shell amplitudes, this last contribution is non-vanishing only when the exchanged particle is on-shell, where it is singular and proportional to δ(u jk − μ 2 jk ). The implications of unitarity for the analytic properties of the 2 → 2 amplitude are well known [47]. In the physical region the partial wave expansion converges and reduces the integral relation given by Eq. (6) to a countable set of algebraic ones. Here s j is the angular momentum of the two-particle system j defined in the IRF jframe, N 2 s j = (2s j + 1)/4π is a normalization constant, f s j (σ j ) is the partial wave amplitude, and P s j ( q j · q j ) is the Legendre polynomial describing the rotation dependence in terms of the cosine of the 2 → 2 scattering angle. The unitarity relation is diagonalized to the partial wave unitarity relation, This equation is automatically satisfied by where the K -matrix is a real function along the unitarity cut. The 3 → 3 amplitude in the physical region can be expanded in partial waves in any of the (12), (13), (23) subsystems. We refer to a subchannel of choice, e.g. (12) as the direct channel and to the others as the cross channels. Since each term in the partial wave expansion is analytic in the angular variables, and therefore in the (13) and (23) invariant masses, singularities in the latter variables can happen only when the series diverges. In contrast to the 2 → 2 case, the unitarity equations for each partial wave would not decouple, and would contain an infinite number of terms. Since in practice one must truncate the series, the amplitude would be regular in the (13) and (23) invariant masses, and the information about the cross channels dynamics would be lost. Instead, we will represent 3 → 3 amplitude in an isobar approximation, where only a finite number of terms in the direct and cross channels are included.

The isobar representation
To be concrete, the partial wave expansion of the connected 3 → 3 amplitude reads where we project the amplitude onto the chosen j and k initial and final channels. Here s j (s k ) is the angular momentum of the initial (final) pair, j ( k ) is the angular momentum between the pair and the spectator, J and M are the total angular momentum of the three particles and its projection, and M J k s k ; j s j is the partial wave amplitude. The angles P j and P k are the orientations of the initial and final pair, which are related to the CMF scattering angle via cos k j = P k · P j . The functions Z J M s contain the rotational dependence of the amplitude A, which are defined as The Z -functions contain all the angular dependence, and they fulfill the orthonormality condition More details are in "Appendix C".
We next discuss the relation between partial wave expansion, isobar representation, and finally the isobar approximation. The partial wave expansion given by Eq. (12) is in principle an exact representation of the amplitude in the physical region of 3 → 3 scattering. However, unlike the analogous expansion in 2 → 2 scattering, the partial wave expansion cannot be used in practice in the 3 → 3 case.
In practice, one needs to restrict the series to a finite number of partial waves. In the physical region of 2 → 2 scattering, the low-energy behavior of the partial waves is determined by barrier factors due to the finite range of interactions. This suppresses the strength of higher partial waves at threshold, provided the latter are regular in the cross channel Mandelstam variables. Cross channel exchanges generate singularities that spoil the convergence of the partial wave series. However, in the 2 → 2 kinematics, these singularities do not overlap with the direct channel physical region. Therefore, the partial wave series can be safely truncated in a finite domain of CMF energies above the two particle threshold. This is not the case, for example, when one of the particles can decay to the other three, and similarly it is never the case for 3 → 3 scattering. If we consider indistinguishable particles, explicit Bose symmetry is lost for the M J k s k ; j s j partial waves, since the partial wave expansion in the initial and final states singles out specific two-body channels. The symmetry is only recovered upon resummation. The isobar representation, in principle, takes care of this problem. One writes the connected 3 → 3 amplitude as a redundant sum of expansions in all the initial and final pairs to make the symmetry explicit. Bose symmetry is thus preserved upon truncation.
As discussed above, one can manage only a finite number of terms in the sums over the subchannel spins. Therefore one reduces the isobar representation to the isobar approximation, by representing the connected 3 → 3 amplitude as a sum over a finite number of isobarspectator amplitudes, as shown in Fig. 6. The truncation is reflected by "max" in the sums. We projected the isobar-spectator amplitudes onto the total angular momentum J of the three particle system. In the following, we refer to A J k s k ; j s j (σ k , s, σ j ) as the partial wave isobar spectator (PWIS) amplitudes. We emphasize Fig. 6 Diagrammatic representation of the isobar approximation amplitude in Eq. (15). The double lines with the black disk represents the isobar amplitude f s j (σ j ), while the gray disk represents the isobarspectator amplitude A k j (p ; p) that, while truncation in s k and s j cannot be avoided in practice, unitarity is diagonal in the total angular momentum. Amplitudes for each J are thus independent and can in principle be resummed. We also stress that the PWIS amplitudes A J k s k ; j s j are not the genuine 3 → 3 partial wave amplitudes M J k s k ; j s j in Eq. (12): where X J k s k ; j s j contains all the cross channel terms which recouple to the direct channel amplitude, The kinematic relations given in "Appendix A" can be used to write the cross channel variables in terms of the direct channel variables. Often in the literature, Bose symmetry is considered as a motivation for Eq. (15). However, this is completely independent: the representation can be applied to the distinguishable particle case (in this case the various A k j (p ; p) contain different physics and have different functional forms), and Bose symmetry can be imposed to the expansion in Eq. (12) without requiring an explicit sum over channels. For example we consider the π + π − π 0 → π + π − π 0 process in the isoscalar vector channel, where the ω is observed. Thinking in isospin basis, where the three pions are indistinguishable, and in the charge basis, where they are distinguishable, leads to the same form of the amplitude, showing that Bose symmetry plays no role in defining the representation.
Isobars parameterize the 2 → 2 dynamics in a given subchannel and angular momentum state. Contrary to the 2 → 2 partial waves, they have only right hand singularities constrained by unitarity. The isobars can have resonant or bound state poles. In the N /D formalism, the isobars can be identified with the denominator function, where the left hand cuts are removed via a dispersive integral [51]. In the following, we will ignore all left hand singularities of the 2 → 2 amplitudes, and identify their partial waves with the isobars. Although we do not need to assume any resonant content for the isobars (e.g. we could use an isobar to describe the π + π + dynamics), it is a popular picture to think of them as a quasi-particle, and to identify the invariant mass and angular momentum of the pair with the isobar mass and spin. Isobars are customarily labeled with the name of the dominant resonance, if any. Isobars can be parameterized as in Eq. (11). For example, the a 1 (1260) decays into three pions dominantly in the ρπ and σ π channels [52]. If one chooses to perform a truncated partial wave expansion of the 3 → 3 amplitude in only the ρπ → a 1 (1260) → ρπ channel, rescattering effects between the ρπ → σ π isobars are ignored. The isobar approximation corrects this by including amplitudes for The approximation is expected to be valid at low values of energy, where a finite number of singularities dominate the amplitude. Moving to higher energies, the left hand cuts controlling the crossed 2 → 4 processes will become relevant, and the behavior of the amplitude will be controlled by analyticity in angular momentum, rather than direct-channel unitarity.
Since the isobar approximation includes the cross channel effects in the summation, the isobar-spectator amplitudes contain only normal threshold singularities determined by unitarity. Therefore, the analytic structure of each isobarspectator amplitude in the energy variables, s, σ j , and σ k , are determined by unitarity.
The problem of convergence in J is more severe. The 3 → 3 amplitude contains an OPE process (see the last diagram in Fig. 5), which can go on-shell in the direct channel, and results in an interaction of infinite range. In this case the cross channel singularities overlap with the physical region and project onto an infinite number of partial waves. The analytic properties of the projected amplitude are highly nontrivial. We discuss them in detail in Sect. 5. However, since the main goal in this and similar studies of three particle scattering is to identify the spectrum, ultimately one needs to deal with amplitudes of well defined total angular momentum J . In other words, these amplitudes diagonalize unitarity, which is the basis for analytic continuation and identification of complex singularities as resonance poles. For this reason, in the following we will not address the problem of convergence in J .

Unitarity relations
It is advantageous to introduce an amputated PWIS amplitude A J k s k ; j s j , in which the isobar amplitudes are factorized, The amputation reduces the number of terms in the isobarspectator unitarity relation by making use of subchannel unitarity in Eq. (10). However, the amputated PWIS amplitudes still have a non-trivial dependence on the subchannel energies due to rescattering effects. As shown in detail in "Appendix C", combining Eqs. (8), (10), (13), (14), (16), and (19) results in the amputated PWIS unitarity relation where C J n s n ; r s r (σ n , s, σ r ) is a purely kinematical recoupling coefficient between different intermediate state isobars, The recoupling coefficients relate two different orientations of three particles in the same frame [44,49,53]. "Appendix C" contains details on the derivation of the recoupling coefficients from the rotational matrices in Eq. (13). The helicity angles and the CMF angle between particles j and k, θ k j , are functions of the invariants (cf. "Appendix A"). The second term contains two integrals over the Dalitz region of the three-particles in the intermediate state, where the physical region is bounded by σ (th) is a function of σ n and gives the physical boundary cos χ n = ±1, e.g. for n = 1 and r = 3, Equation (20) is illustrated in Fig. 7. "Appendix C" contains a sketch of the derivation of the amputated PWIS unitarity relations. The first term of Eq. (20) involves the direct propagation of an isobar in the intermediate state, whereas the second, third, and fourth term involve the exchange of a particle between isobars. The rescattering between isobars modifies the line shape of the isobar amplitudes [54,55]. The final term is the contribution from the OPE process, which gives and additional imaginary part to the amplitude in the physical region. At this stage we have not factored out the threshold factors from partial waves. This is straightforward to implement, however we do not do it here as we consider angular momenta in S-wave in further sections.

The B-matrix parameterization
Motivated by S-matrix theory, we present a parameterization for the PWIS amplitudes that satisfies real axis unitarity given by Eq. (20). In the 2 → 2 case, the K -matrix, is an example of a parameterization satisfying unitarity. For the 3 → 3 case, we present the Bmatrix parameterization for the PWIS amplitudes. The Bmatrix parameterization is a linear integral equation for the amputated PWIS amplitudes that satisfy the unitarity relations Eq. (20): where the B-matrix B J k s k ; j s j contains two terms, The function E J k s k ; j s j is the amputated partial wave OPE amplitude, R J k s k ; j s j is a real function that represents the Fig. 7 Diagrammatic representation for the amputated PWIS unitarity relation in Eq. (20). The black disks in the internal legs represent the isobars, which are amputated from the external legs, see Eq. (19). The cuts across the OPE in the intermediate states yield recoupling coefficients short-distance three-body interactions unconstrained by unitarity, and τ n is the product of the isobar-spectator phase space between and of the isobar amplitude with The parameterization is diagrammatically represented in Fig. 8. The OPE amplitude is defined as where we note that the OPE only contributes to off-diagonal amplitudes, i.e. j = k. In principle, the OPE could contain a regular function of the energy in addition to the pole term, however unitarity only constrains the pole, and we assume all other real functions to be absorbed by R. The amputated partial wave projected OPE amplitude E J k s k ; j s j can be constructed using Eqs. (16) and (19). By construction, R is defined to have no threshold singularities in the threeparticle physical region. The R represents the freedom of short-distance physics for the scattering of three particles, and can be any real function. The function R contains information on three body resonances, and can be modeled with pole terms similar to the K -matrix in the 2 → 2 case, so long as there are no singularities in the physical region of the amplitude. In an effective field theory approach, it represents a low order polynomial of contact interactions. For simplicity, in the following we assume the latter for R. "Appendix D" illustrates how the B-matrix parameterizations satisfies the amputated PWIS unitarity relations. Aspects of its analytical properties are examined in Sect. 5.
The B-matrix parameterization in Eq. (23) differs from Mai et al. [45] in the lower limit of the integral: the latter is derived using Lippman-Schwinger equations with a relativistic potential model, and includes contributions from the unphysical subthreshold region, σ n < σ (th) n . Obviously, both parameterization have the same imaginary part in the physical region, since both satisfy unitarity.
For notational simplicity, let , so that the amplitudes are matrices in the isobar subenergies and angular momenta, which are indicated by the spectator indices. Equation (23) is then a matrix relation with the integrations over intermediate isobars formally represented as matrix multiplications. Recalling that we work with the convention that isobars exists only in the (12) and (32) channels, we write the B-matrix parameterization as the set of coupled equations with the other two amplitudes given by a similar set of equations, The Eq. (28) can be combined into one integral equation for A 13 , where the kernel K 11 is Similarly, Eq. (29) give where the kernel K 33 is given by exchanging the 1 ↔ 3 indices in Eq. (31). Eqs. (30) and (32) can be formally inverted to yield the solutions, Several terms can be identified in the kernels, k j (s), where G is a bubble diagram, H is a box diagram, and the T 's are triangle diagrams, generated by integrals over OPE and contact terms in Eq. (24). Explicitly, Fig. 9 The denominator of the B-matrix parameterization contains four primitive diagrams associated with the rescattering of the B-matrix: a bubble diagram, b and c triangle diagrams, and d box diagram These diagrams occur in the denominators of the amplitudes in Eqs. (33), cf. Fig. 9. They differ to the Feynman diagrams obtained in a perturbative QFT since the integrations are only over the physical region, changing the analytic structure below threshold (see Sect. 5).
The solutions can be interpreted as an infinite series of exchange and bubble diagrams. For example, expanding the solution for A 13 , The first term is the OPE and contact interaction, the second term is a ladder diagram with three exchanges, and various combinations of bubbles and OPE, and so on. The unitarization of bubble diagrams has been considered in quasi-twobody models [56][57][58][59]. In these models it is easy to show how additional cuts appear in the unphysical sheets due to the isobar decay.
Three-body resonances manifest as poles in the complex s-plane of the scattering amplitude. Rearranging the constituents of the kernel relates the two denominators Thus, we can write the full 3 → 3 amplitude in terms of a single Fredholm determinant. The determinants are independent of the external isobar energies, and the intermediate integrations will modify the phase space factors to incorporate rescattering effects. Resonance poles can be determined by solving The B-matrix solutions are real-boundary values of analytic functions in the complex s-plane. The physical amplitudes are defined by s → s + i , σ j → σ j + i , and σ k → σ k + i . Aspects of its analytic properties are discussed in the following section.
Relation to the finite volume formalism In finite volume studies for lattice QCD, substantial progress has been made to understand the connection between discrete energy levels and properties of hadron scattering amplitudes [29][30][31][32][33][34][35][36][37][38][39]. In the case of 2 → 2 scattering the two-particle finite volume spectrum constrains the values of the infinite volume partial wave amplitudes via the Lüscher quantization condition [60]. The multi-variable nature of 3 → 3 scattering amplitudes makes the derivation of the finite volume quantization condition much more complicated and different groups have approached the problem from a different angle. For example, in Refs. [29][30][31][32][33] the authors introduce amplitudes labeled by subchannel spins and the spectator 3-momenta. Furthermore, ladder diagrams generated by OPE are considered independently from other interactions. This implies that partial wave projection to total spin, which is necessary if one is interested in extracting properties of three-body resonances, would be performed after resummation of the OPE ladder. On the other hand, in Refs. [34,35], the quantization conditions are derived starting from a set of amplitudes projected onto the total and subchannel spins from the start, i.e. before OPE resummation, in a spirit close to our work. In Refs. [36][37][38] the quantization conditions are derived in a nonrelativistic EFT framework, and the direct comparison with our S-matrix approach is more complicated. It is more interesting to discuss the differences with Refs. [29][30][31][32]. Since we do not aim to address the subtleties of the finite volume here, we compare with the infinite volume equations derived in there on the basis of the finite volume formalism. For simplicity we ignore coupling to the 2-body channel. In Refs. [29][30][31][32], the 3 → 3 connected amplitude is denoted by M 33 ( k, k ) (see Eq. (112) in Ref. [32]). It contains the resummed OPE ladder and the amputated amplitude T 33 ( k, k ) that is generated by the kernel K df,33 ( k, k ), which is analogous to our driving term R J k s k ; j s j (σ k , s, σ j ).
Both the OPE ladder and the amplitude T 33 ( k, k ) are solutions of linear integral equations (see Eqs. (87) and (106) in Ref. [32]), which are analogous to our Eq. (23).
To further illustrate the connection between our amplitudes and those of Refs. [29,32] we shall consider the case of three identical particles in S-waves. The phase space ρ 2 (Eq. (7)) of the two particle subsystem gains a factor 1/2! to account for their identical nature. The resulting unitarity relations have a similar form as in Eq. (20). In matrix notation, one finds where the matrices are in the σ , σ space with f and ρ 3 diagonal matrices. The S-wave projection of the OPE is given by the symmetric matrix E, and is found by the inverse relation of Eq. (16) on the OPE amplitude in Eq. (27). It is straightforward to show that the B-matrix parameterization (cf. Eq. (23)), satisfies the unitarity relation, Eq. (38). Moreover, after simple manipulations, Eq. (39) can be rewritten as where the D amplitude is the ladder sum of OPE, given by and L ≡ f + Dρ 3 . Finally we introduce the amplitude T satisfying and obtain an expression closely resembling that in Refs. [29,32], The difference between Eq. (43) and the corresponding expression for M (u,u) 33 in Refs. [29,32] is in the definition of L. In our notation, the L of Refs. [29,32] contains an additional 1/3 constant, and the f and D matrices contain an extra factor of ρ 2 .
Although these analogies should be verified with care, two main differences appear. One is in the treatment of the OPE dynamics, which in Ref. [32] is resummed before projection onto the total spin and in our case the projection is done first. It is likely that these approaches will ultimately prove to be equivalent, since in practical applications only a finite number of partial waves in total spin or spectator momentum components can be kept. The other difference is in the L function, which could possibly be related to a discrepancy in the definitions of R and K df,33 . It would be interesting to see if our approach and the corresponding equation of Refs. [29,32] provide the same result, and to determine the origin of the difference.
One can also consider our formalism in the finite volume by relating integrals over the isobar invariant mass to discrete sums. This approach would lead to a quantization condition, similarly to what was shown by Ref. [34]. It remains to be seen if a quantization condition derived in this manner is identical to that of Ref. [29]. This is an active area of research, however, outside the scope of this work, and we leave it for future studies.

Aspects of analytic properties
In this section, we examine the singularities of the OPE amplitude and the triangle diagram from the B-matrix parameterization. We numerically evaluate an amplitude where all external particles have unit mass (m 1 = m 2 = m 3 = 1) and coupling. In these studies, the units are arbitrary. For simplicity, we consider S-waves only, i.e. J ( s ) k ( s) j = 0(00) k (00) j . Generalizing to nonzero angular momenta does not change the analytic properties.

One particle exchange
As seen in Eqs. (34), the building block for the B-matrix kernels is the OPE amplitude. Projecting Eq. (27) using Eq. (16) gives the S-wave OPE amplitude, where z k j is given as where λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + bc + ca) is the Källén triangle function. Eq. (45). We investigate the OPE as a function of s for fixed real σ j and σ k . The imaginary part of the OPE is which is given by the unitarity relations in Eq. (20). The OPE has four branch points in s, one at zero, one at infinity, and two which we label s k j , called the real particle exchange (RPE) cut. The VPE cut is associated with the exchange of a virtual particle, generating long-range forces. Historically, the RPE cut is associated with the exchange of a real particle between isobars, i.e. when it is kinematically allowed for an isobar to decay. This corresponds to when the RPE branch points lie on the real axis above the isobar-spectator threshold. If the isobar invariant masses are below the decay threshold, then the RPE branch points move in the complex plane below the isobar-spectator threshold. For convenience, however, we will always call this the RPE cut, and emphasize that a real particle exchange occurs only if it is kinematically accessible. Note that although the value of the isobar mass dictates the physics of the OPE, the OPE is blind to the decay products of the isobar and the physical threshold in s is max{( √ σ j + m j ) 2 , ( σ k + m k ) 2 }. We can understand the analytic structure of the OPE by writing a dispersive representation in s. Equation (46) is nonzero in two regions, leading to the relation where V is the contour over the VPE cut and R is the contour over the RPE cut. The integrand has four branch points associated with the thresholds and pseudo-thresholds of the initial and final momenta. We choose to orient the branch cuts such that the lowest branch point (min{( σ k − m k ) 2 , ( √ σ j − m j ) 2 }) has a cut running to −∞, the highest branch point (max{( σ k + m k ) 2 , ( √ σ j + m j ) 2 }) has a cut running to +∞, and the other two branch points have a branch cut joining them. The contour V is always taken above the real axis, whereas the contour R depends on the external masses. The physical amplitude is defined as the boundary value when s → s + i , below the RPE cut.
For fixed σ j > σ (th) j , the RPE cut can be categorized by different regions in σ k . In Fig. 10, we illustrate how the analytic structure of the integrand and the integration contours change in these regions. Assuming a small imaginary part μ 2 jk → μ 2 jk − i , the RPE branch points have a finite imaginary part for σ k < σ (th) k , with opposite signs. When σ k > σ (th) k , the branch points are infinitesimally close to the real axis. In the physical region, s (−) k j has inflection points at two locations of σ k : which follow from ds We can therefore classify the regions according to when the RPE branch points are both in the upper-half plane or when they approach the real axis. k j is above the real axis. The RPE cut connects these two points by crossing the real axis below the threshold ( σ k + m k ) 2 . Real particle exchange in this case has consequences when considering the OPE processes embedded in the triangle diagram, which is discussed in the next section. For k = 1 and j = 3, Fig. 10a shows the RPE and VPE contours, R and V , respectively. Note that in the logarithmic representation, the RPE cut is circular, whereas in the dispersive representation, one can define the cut in any chosen manner as long as singularities are not crossed. Fig. 11. When σ k decreases below the inversion point σ k j without crossing the real axis. This is the typical case when considering the exchange of a real particle, illustrated in Fig. 10b. Note that when σ k = σ j for equal masses m j = m k , then the integrand branch points merge into pole singularities.
k , see region (c) in Fig. 11. The RPE branch points again wrap around the real axis, cf.
1 , and d σ 1 ≤ σ (th) 1 . Real particle exchange cannot occur in case (d). In the logarithmic representation Eq. (44), the RPE cut is circular k , see region (d) in Fig. 11. The branch points s (±) k j move deep into the complex plane, as shown in Fig. 10d. In this region, the isobar cannot decay, and therefore it is unphysical for the 3 → 3 elastic scattering. Real particle exchange cannot occur, leaving only the virtual contributions.
If we evaluate the OPE along the real s-axis in regions (a) or (c), we find that the real part of the OPE has a jump due to crossing the RPE cut. In the logarithmic representation, this crossing occurs when z k j = 0, that is, when These cases are illustrated in Fig. 10 for spectators k = 1 and j = 3. We plot the OPE amplitudes, Eq. (44), as a function of s for fixed σ 3 and σ 1 in Fig. 12. Figure 12a shows the OPE computed at σ 1 in region (a). At this energy, the s (−) k j is below the real axis, and the RPE cut wraps around the real axis, passing below the threshold ( σ 1 + m 1 ) 2 . The jump in the real part at s (0) 13 is due to crossing the RPE cut. Figure 12b, is evaluated at σ 1 in region (b), where both branch points are above the real axis. Here, we illustrate that as σ 1 decreases, the width of the imaginary part decreases and the peak increases. The narrowing imaginary region physically represents that less phase space is available for real particle propagation in the intermediate state. Figure 12c is computed for σ 1 in region (c), right above the two-particle threshold. There is a jump in the real part at s (0) 13 from crossing the RPE cut. The final case is illustrated in Fig. 12d, where the OPE computed in the unphysical region (d). There is an imaginary part due to the VPE cut only, as it is kinematically inaccessible for the exchange of a real particle. The jump in the real part at s The physical region is taken as the region approaching the real axis, below the RPE cut.
To summarize, the analytic structure of the OPE is given by two branch cuts, the VPE and RPE cuts. The VPE cut is present for −∞ < s ≤ 0, and is associated with the exchange of an off-shell particle. For physical isobars, the RPE cut is in the physical region. We have shown different scenarios, identified by the isobar masses, in which the RPE branch points can approach the physical region, which impact the structure of the B-matrix kernels.

Triangle diagrams
To understand resonance poles of 3 → 3 systems, the analytic structure of the B-matrix parameterization, Eq. (33) must be understood in the complex s-plane. This means understanding the properties of the B-matrix kernels. Here, we investigate the triangle diagram, and leave the box diagram for future studies. Let us work with the triangle T B ≡ T (2) 11 introduced in Eq. (34c), where all angular momenta are in S-wave. For convenience, let R = 1, thus the amplitudes are independent of σ 1 , and given by where , and the dependence of T B on σ 1 has been understood. To ensure the correct analytic properties of the isobar amplitude, we introduce its dispersive representation giving the form for T B We see the σ 3 -integral does not depend on f 3 ( σ ), so for simplicity we take the narrow width limit Im where M is the mass of the isobar. The narrow width limit shifts the unitarity cut in the triangle diagram to begin at the threshold s = (M + m 3 ) 2 . However, for a general isobar shape, Eq. (53) can be used to sum over its distribution, recovering the correct unitarity branch cut starting at s th . Therefore, the triangle diagram has the form Figure 13 shows the triangle diagram in consideration. The B-matrix triangle contains singularities on the physical sheet. These are due to the s-singularities in ρ 3 and E S 13 , and to endpoint singularities when the integration limits hit the σ 3singularities of the integrand. The upper integration limit gives a branch cut for s < 0. Since ρ 3 (s, σ 3   On the real s-axis, the OPE has a discontinuity in the real part when z 31 = 0, i.e. at σ 3 = σ (th) 3 . This discontinuity from crossing the RPE cut is present in T B . Two more singularities occur when σ 3 hits the two inflection points σ (a) 3 and σ (b) 3 , which are defined in Eq. (49). The OPE pinches the real axis at σ 3 = σ (a) 3 , and generates a singularity in T B at the initial state threshold s = ( √ σ 1 + m 1 ) 2 . This can be understood by realizing that the OPE branch points, Eq. (47), can alternatively be written in terms of σ 3 as a function of s for fixed σ 1 . The branch points are then σ and σ (−) 3 lies infinitesimally below the real axis in the physical region. Figure 14 shows the motion of σ (±) 3 in the complex σ 3 -plane as a function of s for fixed σ 1 . At the three particle threshold, s = (m 1 + m 2 + m 3 ) 2 , the branch points have finite imaginary part and are on opposite sides in the σ 3 -plane. As s approaches the initial state threshold s = ( √ σ 1 + m 1 ) 2 , the σ (±) 3 branch points pinch the real axis at σ 3 = σ (a) 3 . Since the T B integration is on the real axis starting from σ (th) 3 , the integration path is pinched, causing a singularity in T B at s = ( 3 , the branch point migrates back below the real axis at a 3 , this effect generates the triangle singularity [50,61,62]. The triangle singularity has been studied as a possible mechanism to explain anomalous structures observed in heavy flavor experiments [19,[63][64][65]. The peak of the triangle singularity coincides with the s (−) 31 branch point, i.e. s tr = s (−) 31 . Aside for the unitarity branch cut starting at s = (M + m 3 ) 2 and the triangle singularity, these additional singularities in the physical s-plane are not allowed by analyticity. The extra singularities are moved to the second sheet when we consider the integration over the isobar shape, cf. Eq. (53), leaving only the unitarity cut starting at s = (m 1 +m 2 +m 3 ) 2 and the triangle singularity.
We compare the structure of Eq. (54) with that of a Feynman diagram triangle in a QFT (see "Appendix E" for a review of the Feynman triangle), which can be written as where T is the path from the threshold (M +m 3 ) 2 to ∞, and the S-wave amplitudes are normalized according to Eq. (16). Figure 15 shows the real and imaginary parts, respectively, of the two triangle diagrams T F and T B , below the region of the triangle singularity. Notice that the Feynman triangle has only a normal threshold singularity at s = (M + m 3 ) 2 , and is smooth everywhere else. The imaginary parts of both triangles are identical above threshold, as required by uni-tarity. The B-matrix triangle has noticeable kinks in both the real and imaginary parts below threshold, corresponding to the singularities discussed above. The black dashed lines indicate the location of the singularities. Starting from low energy, the first additional singularity is the s = 0 singularity from the phase space. The next two singularities occur at 3 + m 3 ) 2 , which are from the phase space evaluated at the lower integration limit. The real part contains a singularity from evaluating the OPE across the RPE cut. Note that the imaginary part does not contain this jump, consistent with the OPE description in the previous section. The next singularity occurs at the initial state threshold s = ( √ σ 1 + m 1 ) 2 , which is due to the pinching of the σ 3 contours by the OPE branch points. Finally the normal threshold at s = (M + m 3 ) 2 . Figure 16 shows T B and the OPE in the region where the triangle singularity develops. The line shape shows the production threshold at (M +m 3 ) 2 , and the peak at s = s tr . The OPE branch point s (−) 31 clearly coincides with the triangle peak. Figure 17 shows the T B and T F as a function of s at fixed σ 1 and varying M 2 in the region below and above M 2 = σ (b) 3 . Figure 17a shows the real parts, and Fig. 17b shows the imaginary parts. At  Figure 17c shows the RPE cut of the OPE in the s-plane at the corresponding values for the triangle amplitude. We also compare the B-matrix triangle with the analogous one from Mai et al. [45], that we denote as T M , where we take their contact term equal to unity, and the lower integration limit in their model accounts for the physics in the unphysical region. As σ 3 → −∞, the OPE amplitude goes like 1/σ 3 , while the phase space grows as σ 3 , thus the integrand goes like 1/σ 3 and the function is logarithmically divergent. Numerically, we choose to cut off the integral at some large value, e.g. −200 to investigate the behavior. For T M , all lower limit endpoint singularities in s from the phase space and OPE are moved toward −∞. The s = 0 pole from the phase space persist, and the normal threshold singularity at s = (M + m 3 ) 2 is present since it is from the upper limit.
The pinch singularity at s = ( √ σ 1 + m 1 ) 2 is also present, as well as the pinch singularity at s = ( √ σ 1 − m 1 ) 2 , cf. Fig. 14. The second pinch singularity occurs when the integration over σ 3 hits σ 3 is a third inflection point in the unphysical region, occurring at s is at the threshold s = ( √ σ 1 − m 1 ) 2 (when Im s (−) k j = 0). This pinch singularity is absent in the B-matrix triangle, as the integral is only over the physical region. Figure 15 compares the line shapes of all three triangles, T B , T F , and T M . Although T M has a logarithmic divergence, we fix the lower integration limit to −200. We see how the line shape below threshold smooths out except at the remaining singularities, shown with the black dashed lines. The red dashed line indicates the second pinch singularity in T M .
The Feynman triangle can be recovered from T M with the method discussed by Aitchison and Pasquier [42], where the isobar approximation for 1 → 3 decays was studied. Using their inversion technique, it was found that the Feynman triangle can be written as a dispersive integral over the isobar invariant mass as in Ref. [45], plus additional terms. The latter are real in the physical region, but cure the below threshold singularities shown in the B-matrix. The additional terms also cancel the logarithmic divergence, leaving a finite amplitude. The addition of terms which cancel the unphysical singularities can be traced to what is chosen for R. The R-matrix has the freedom to absorb these differences as they modify the structure below threshold, leaving the physical region unchanged.

Removal of unphysical singularities
As shown, the B-matrix parameterization contains additional singularities which do not match the expected analytic behavior of the amplitudes. This happens in both our formulation and Mai et al. [45]. The formalism only considers unitarity, which constrains the singularities in the physical region. Imposing additional constraints from analyticity would in principle remove these extraneous singularities, as in the 2 → 2 case. One possible venue for improving that is to substitute the B-matrix kernels, Eq. (34), with the Feynman one. This is in the same spirit of the Chew-Mandelstam phase space in the 2 → 2 parameterizations, which removes the unphysical singularities of the phase space. Although the kernels will now have the proper analytic structure (no physical sheet singularities except for the unitarity cut), the resulting amplitude will still contain singularities from iterating the kernel. Consider the solution for A 33 in Eq. (33b), where the kernel is replaced by the Feynman one, Now expand the solution Eq. (33b) in an infinite series, The first term is the kernel, composed of Feynman diagrams which have the correct analytic properties. Let the kernel consist only of the triangle diagram, then the second term is two Feynman triangles joined with a τ -function. The equivalent Feynman diagram would have two exchanges integrated over the four-momenta, which is not equivalent to what is shown in Eq. (60) due to the τ -function. This diagram, as well as the higher-order ones, contain non-analyticities in a similar manner to what was shown for the triangle diagram. The unintegrated singularities from the phase space are always present. Therefore the simple kernel substitution does not produce the correct analytic behavior in the B-matrix solution. However, it can still be advantageous, as it corrects some of the unphysical singularities in the present B-matrix solution.
The remaining singularities should disappear if one was to solve the proper Bethe-Salpeter equations of the underlying QFT. The B-matrix parameterization is indeed reminiscent of that for 2 → 2 scattering. We examine some differences between these formalisms. The B-matrix parameterization is a covariant integral equation for the on-shell isobar-spectator amplitudes. It satisfies unitarity relations and does not have additional imposed constraints from analyticity. Thus, for complex energies on the physical Riemann sheet, the Bmatrix parameterization contains the unitarity cut, and has additional s-singularities from the τ and OPE.
The Bethe-Salpeter equation is a covariant integral equations that incorporate an infinite number of exchanges for any given QFT [48]. Solving it amounts to summation of exchange diagrams, similarly to the B-matrix. The resulting amplitudes are analytic functions in the complex s-plane, as the QFT amplitudes inherently obey analyticity constraints. The physical sheet thus has only the allowed singularities, such as the unitarity cut and possible bound state poles. Lippmann-Scwhinger equations are nonrelativistic equations for the scattering amplitude in a given potential model. The B-matrix has similarities to the Lippmann-Schwinger equations in that both involve in a three-dimensional integral over the momenta [45]. In this work however, we focus on the physical region, and truncate the isobar mass integration appropriately. Conversely, the Bethe-Salpeter equation contains integrations over four-momenta, which results in integrating over the off-shell behavior of the amplitude. Introducing dispersive integrations in the B-matrix amounts to the same procedure, and would remove the unphysical singularities.

Conclusions
In summary, we have discussed the phenomenological description of 3 → 3 elastic scattering of spinless particles. The 3 → 3 amplitude was described in the isobar representation. We constructed the unitarity relations for the isobarspectator amplitudes for general partial wave quantum numbers. For a practical use, the infinite sums are truncated, leading to the standard isobar approximation. We parameterize the isobar-spectator partial wave amplitudes with the Bmatrix formalism, which automatically satisfies the unitarity. The B-matrix parameterization explicitly includes the one pion exchange as a long-range contribution required by unitarity. The short-range part is not constrained by unitarity, and it can be incorporated by a specific (model-dependent) choice of the parameterization. This gives to the framework enough freedom to incorporate QCD resonances. The approach here differs from Mai et al. [45] in that the 2 → 2 amplitudes required as input are only needed to be known in the physical energy regions. The singularities of the OPE directly impact the analytic structure of the B-matrix kernels, and are discussed explicitly for the triangle-like diagram. The singularities in the unphysical region of our solution differ from the Mai et al. ones, and from the Feynman diagram triangle. This results in a different value for the real part of the amplitudes in the physical region. Further studies are needed to understand how to remove unexpected singularities from the B-matrix. We also compare our formalism to the most recent ones discussed in the literature to extract three-body scattering amplitudes from lattice QCD. In particular, the main difference with Refs. [29][30][31][32][33] consists in the order of how the partial wave expansion and the one particle exchange ladder summation is performed. It remains to be seen whether the two operations commute, and whether the resulting amplitudes coincide.
Future studies will investigate the continuation to the unphysical energy sheets. This venue will allow us to constrain the role of one particle exchange in generating resonant structures, as it is assumed in some molecular models for the XY Z states.
where μ jk is the mass of the particle that is neither j nor k.
The remaining variables needed to completely describe the 3 → 3 process are found by examining the IRFs. The initial IRF j and final IRF k are defined when P j = 0 and P k = 0, respectively. We use the convention that initial and final state variables are evaluated in their own respective IRF. The momentum of the first particle in the initial pair is denoted as q j in the IRF j . Similarly, the first particle in the final pair is q k in the IRF k . For example, for the final spectator 3 in the IRF 3 , q 3 is the final momentum of particle 1, and in the IRF 1 of spectator 1, q 1 is the initial momentum of particle 3. In terms of invariants, these momenta are The spectator momenta in these frames are for the final state and for the initial state. The helicity angles of the first particle in the IRFs are given by χ j and χ k , for the initial and final states, respectively. The helicity angles are defined w.r.t. the opposite line-of-flight of the spectator. The azimuthal angles for the initial and final state are γ j and γ k , respectively. The azimuthal angles are defined as the angle between the plane of the two particles in the CMF, and the IRFs, cf. Fig. 3. Note that the azimuthal angles γ j and γ k are invariant with respect to the Lorentz boost between the IRFs and CMF, so γ j = γ j and γ k = γ k . The invariant masses of the other two pairs in their respective frames are related to the helicity angles. For example, in the IRF 3 , (A7)

Appendix B: Derivation of 3 → 3 unitarity relations
In this appendix we derive the general elastic unitarity relations for the 3 → 3 elastic scattering of distinguishable spinless particles [40,41,50]. For convenience, in this section we adopt the notation that the normalization of a single particle state is p k |p j = (2π) 3 2ω j δ (3) The S-matrix is a unitary operator, S † S = 1, which implies that T − T † = i T † T , where S = 1 + i T . We consider the system in an energy range above the three particle threshold, but below the first inelastic threshold, s th ≤ s < s inel . Taking matrix elements of this operator between initial and final states |p and |p , and inserting the completeness relation 1 = d p 1 d p 2 d p 3 |p p |, gives the unitarity relation The matrix elements p | T † |p are equal to p | T |p * by the property of Hermitian analyticity [50,66]. Thus the left hand side of Eq. (B1) gives the imaginary part of the matrix element, The right hand side of Eq. (B1) is evaluated by substituting Eq. (B2) and expanding the product into four terms, The fourth term contains two cases, one where j = k, and one where j = k, so we split the sum into the two distinct terms We can write δ( p j − p j ) δ( p j − p j ) = δ( p j − p j ) δ( p j − p j ) in the first term in Eq. (B5), thus we can identify the disconnected unitarity relation as being proportional to the spectator singularity δ( p j − p j ), and the connected unitarity relation The momenta with j 1 and j 2 in Eq. (B6) identify the first and second particle in the pair. Substituting Eqs. (2) and (3) into Eqs. (B6) and (B7), and evaluating the phase space integrals yield the unitarity relations Eqs. (6) and (8).

Appendix C: Derivation of PWIS unitarity relations
Using the assumption of the isobar model Eq. (15), we derive a set of unitarity relations for the amputated PWIS amplitudes. Inserting Eq. (15) into the unitarity relations Eq. (8) leads to a unitarity relation for the A k j isobar-spectator amplitude, where we wrote the three-body phase space factor in the first two terms, where P n is the orientation of the isobar associated with spectator n in the intermediate state, and q n is the orientation of the first particle in the intermediate isobar in its rest frame.  Figure 18 shows the diagonal and off-diagonal topologies in the intermediate state.
The partial wave expansion Eq. (16) can be derived by considering the expansion in three steps. The first step is to perform the expansion of the isobars into definite spin and helicity, where λ j , λ k are defined along the direction of the isobar in the CMF. The expansion removes the q k and q j dependence in the helicity amplitude A s k λ k ;s j λ j . Second, the helicity amplitude can be expanded in partial waves, where N 2 J = (2J + 1)/4π . Finally, since parity is not a good quantum number in the helicity basis, we convert the helicity partial wave amplitudes into L S partial wave amplitudes, where C J sλ = √ (2 + 1)/(2J + 1) J λ| 0sλ , and has the completeness relation Combining Eqs. (C3), (C4), and (C5) yields the partial wave expansion Eq. (16). We apply the expansion to Eq. (C1) to obtain the PWIS unitarity relations. The diagonal terms are most directly evaluated using the orthonormality condition, Eq. (14). Since p k = p k in the third term, and p j = p j in the fourth term, then the intermediate isobar orientation is identical to that of the final and initial state isobar, respectively. The integrations over q k and q j can be performed by writing Eq. (9) using the spherical harmonic addition theorem. The off-diagonal terms are more challenging, as they involve two different angles in the intermediate state, thus the rotational functions will not directly integrate out. We can make use of the group properties of rotations to simplify the intermediate rotational functions to a recoupling coefficient. The off-diagonal terms on the right-hand side of Eq. (C1) under the expansion Eq. (16) will contain terms of the form where n = r , and we combined the terms with γ n and the orientation of the isobar to the set of angles R n = (α n , β n , γ n ), where α n is the azimuthal angle of the isobar and β n is the polar angle, w.r.t. some fixed coordinate system. Note that since we boost along the direction of the isobar to go between CMF and IRFs, γ n = γ n . The angles R are the Euler angles describing the orientation of the three particles in their CMF. Since these two sets of angles describe the same configuration of three particles, with the only difference being which particle is the spectator, the angles R n and R r must be related by a rotation. Each set of angles can be found by rotating from some initial standard configuration. We define the standard configuration such that the three particle system lies in the xz-plane, where the spectator momenta is along the negative z-axis, cf. Fig. 19. Then, the difference in the Euler angles is a rotation about the y-axis, R r = R n r nr , where r nr is the rotation relating the two standard configurations [44,49]. Here, the rotation is about the y-axis, r nr = R y (θ nr ), where θ nr is given in Eq. (A3). Thus, the rotation is a function of the invariant masses of the isobars, θ nr = θ nr (σ n , s, σ r ). Note that the inverse rotation is given by r rn = r −1 nr .
(a) (b) Fig. 19 The standard configurations considering a particle 1 as the spectator and b particle 3 as the spectator Therefore, we can relate the two Wigner D-matrices using the group addition property where d R n = dα n d cos β n dγ n . Since the angle χ can be written in terms of the invariant masses of the isobars, it is advantageous to write the phase space in terms of the Dalitz variables, where we used The last piece needed is the partial wave projection of the OPE term. To evaluate the partial wave projection, we write the delta-function as where z k j is defined in Eq. (45). Then, the completeness relation for the delta-function allows us to write Eq. (C14) as where λ and λ are arbitary, and thus we may choose them to align with λ j and λ k , respectively. Then, d Finally, the 2 → 2 amplitudes are written using Eq. (9) whereq j is the momentum of the first particle in the final state in the IRF j , andq k is the momentum of the first particle in the initial state in the IRF k . Since p k = p k and p j = p j , then the azimuthal angles are identical, γ k = γ j . The helicity angles of the first particle in the opposite frames are defined as χ j and χ k , cf. Fig. 20. However, we can easily identify that χ k = χ k and χ j = χ k since the intermediate spectator is aligned for the OPE and the IRFs merely differ by the rotation about z.
The imaginary part of the kernel is Grouping common terms in Im τ n and Im B k j , and identifying the forms of the amplitudes from Sect. 4 Substituting for the imaginary parts of τ n and B k j gives the PWIS unitarity relation for A 13 . The unitarity relations for the other amplitudes can be found in a similar manner.

Appendix E: The Feynman triangle diagram
For reference, we state the basic formulae for computing the Feynman triangle diagram, cf. Ref. [48]. The perturbative Feynman diagram has the form shown in Fig. 21, where the denominator is the product of internal propagators, Using the Feynman parameterization and standard loop integration techniques, the Feynman diagram has the form