Solving relativistic three-body integral equations in the presence of bound states

We present a systematically improvable method for numerically solving relativistic three-body integral equations for the partial-wave projected amplitudes. The method consists of a discretization procedure in momentum space, which approximates the continuum problem with a matrix equation. It is solved for different matrix sizes, and in the end, an extrapolation is employed to restore the continuum limit. Our technique is tested by solving a three-body problem of scalar particles with an $S$ wave two-body bound state. We discuss two methods of incorporating the pole contribution in the integral equations, both of them leading to agreement with previous results obtained using finite-volume spectra of the same theory. We provide an analytic and numerical estimate of the systematic errors. Although we focus on kinematics below the three-particle threshold, we provide numerical evidence that the methods presented allow for determination of amplitude above this threshold as well.


I. INTRODUCTION
Several outstanding problems in modern-day hadronic, particle, and nuclear physics require a relativistic description of the dynamics of multi-hadron systems. Many resonances, which challenge our understanding of the strong interaction, are observed experimentally in reactions involving final states composed of three particles or more. One example is the recently observed tetraquark candidate X(2900) found in the B + → D + D − K + decay [1,2]. Due to the complexity of these reactions, it is rarely evident if these are indeed genuine states of quantum chromodynamics (QCD), or merely kinematic enhancements [3,4]. Similarly, three-body decays play a significant role in modern-day tests of the fundamental symmetries of the Standard Model and searches of its extensions. A prominent example is the measurement of the enhanced CP violations in B ± decays to three light mesons [5], where the large CP asymmetries can result from the presence of a rich resonant structure in the three-body final state. Lastly, it is well known that the three-nucleon forces are indispensable in the effective description of light nuclei and their properties [see Refs. [6][7][8] for some key examples]. However, the exact form of their contribution, within the context of QCD, is still undetermined, see Ref. [9].
To resolve these, and many other problems, a coordinated effort has been initiated to obtain two-and three-hadron dynamics from QCD using lattice QCD. 1 Although the scattering amplitudes are not accessible directly in the finite volume computations, one can obtain them from the finite-volume spectra computed with lattice QCD via appropriate non-perturbative mappings called quantization conditions [12,13]. This technique has proven successful in the twohadron sector [14][15][16][17][18][19][20][21][22][23][24][25][26][26][27][28][29][30][31][32][33][34], including systems where multiple open channels are kinematically accessible [35][36][37][38][39][40][41][42][43][44]. Extensions of this methodology to the three-particle sector have been formally developed in recent years, focusing on three identical scalar particles. Two approaches have been followed to address the determination of relativistic three-particle scattering amplitude from lattice QCD. The first is the relativistic field theory (RFT) approach, which derives the quantization condition by summing on-shell projected generalized Feynman diagrams to all-orders [45][46][47][48][49][50][51][52][53][54]. An alternative method based on S matrix unitarity, called the finite volume unitarity (FVU) approach, constructs on-shell scattering equations from the unitarity relation for the amplitude [55][56][57][58], then postulates a finite volume analog which can be used to derive a quantization condition [59][60][61]. These methods have been shown to result = < l a t e x i t s h a 1 _ b a s e 6 4 = " D m N M k 2 Y X Y y T v + z N 0 r W X Y w O / Y u L g = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 I t Q 9 O K x B f s B b S i b 7 a R d u 9 m E 3 Y 1 Q Q n + B F w + K e P U n e f P f u G 1 z 0 N Y H A 4 / 3 Z p i Z F y S C a + O 6 3 0 5 h b X 1 j c 6 u 4 X d r Z 3 d s / K B 8 e t X S c K o Z N F o t Y d Q K q U X C J T c O N w E 6 i k E a B w H Y w v p v 5 7 S d U m s f y w U w S 9 C M 6 l D z k j B o r N W 7 6 5 Y p b d e c g q 8 T L S Q V y 1 P v l r 9 4 g Z m m E 0 j B B t e 5 6 b m L 8 j C r D m c B p q Z d q T C g b 0 y F 2 L Z U 0 Q u 1 n 8 0 O n 5 M w q A x L G y p Y 0 Z K 7 + n s h o p P U k C m x n R M 1 I L 3 s z 8 T + v m 5 r w 2 s + 4 T F K D k i 0 W h a k g J i a z r 8 m A K 2 R G T C y h T H F 7 K 2 E j q i g z N p u S D c F b f n m V t C 6 q n l v 1 G p e V 2 m 0 e R x F O 4 B T O w Y M r q M E 9 1 K E J D B C e 4 R X e n E f n x X l 3 P h a t B S e f O Y Y / c D 5 / A I z N j M E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D m N M k 2 Y X Y y T v + z N 0 r W X Y w O / Y u L g = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 I t Q 9 O K x B f s B b S i b 7 a R d u 9 m E 3 Y 1 Q Q n + B F w + K e P U n e f P f u G 1 z 0 N Y H A 4 / 3 Z p i Z F y S C a + O 6 3 0 5 h b X 1 j c 6 u 4 X d r Z 3 d s / K B 8 e t X S c K o Z N F o t Y d Q K q U X C J T c O N w E 6 i k E a B w H Y w v p v 5 7 S d U m s f y w U w S 9 C M 6 l D z k j B o r N W 7 6 5 Y p b d e c g q 8 T L S Q V y 1 P v l r 9 4 g Z m m E 0 j B B t e 5 6 b m L 8 j C r D m c B p q Z d q T C g b 0 y F 2 L Z U 0 Q u 1 n 8 0 O n 5 M w q A x L G y p Y 0 Z K 7 + n s h o p P U k C m x n R M 1 I L 3 s z 8 T + v m 5 r w 2 s + 4 T F K D k i 0 W h a k g J i a z r 8 m A K 2 R G T C y h T H F 7 K 2 E j q i g z N p u S D c F b f n m V t C 6 q n l v 1 G p e V 2 m 0 e R x F O 4 B T O w Y M r q M E 9 1 K E J D B C e 4 R X e n E f n x X l 3 P h a t B S e f O Y Y / c D 5 / A I z N j M E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D m N M k 2 Y X Y y T v + z N 0 r W X Y w O / Y u L g = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 I t Q 9 O K x B f s B b S i b 7 a R d u 9 m E 3 Y 1 Q Q n + B F w + K e P U n e f P f u G 1 z 0 N Y H A 4 / 3 Z p i Z F y S C a + O 6 3 0 5 h b X 1 j c 6 u 4 X d r Z 3 d s / K B 8 e t X S c K o Z N F o t Y d Q K q U X C J T c O N w E 6 i k E a B w H Y w v p v 5 7 S d U m s f y w U w S 9 C M 6 l D z k j B o r N W 7 6 5 Y p b d e c g q 8 T L S Q V y 1 P v l r 9 4 g Z m m E 0 j B B t e 5 6 b m L 8 j C r D m c B p q Z d q T C g b 0 y F 2 L Z U 0 Q u 1 n 8 0 O n 5 M w q A x L G y p Y 0 Z K 7 + n s h o p P U k C m x n R M 1 I L 3 s z 8 T + v m 5 r w 2 s + 4 T F K D k i 0 W h a k g J i a z r 8 m A K 2 R G T C y h T H F 7 K 2 E j q i g z N p u S D c F b f n m V t C 6 q n l v 1 G p e V 2 m 0 e R x F O 4 B T O w Y M r q M E 9 1 K E J D B C e 4 R X e n E f n x X l 3 P h a t B S e f O Y Y / c D 5 / A I z N j M E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D m N M k 2 Y X Y y T v + z N 0 r W X Y w O / Y u L g = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 I t Q 9 O K x B f s B b S i b 7 a R d u 9 m E 3 Y 1 Q Q n + B F w + K e P U n e f P f u G 1 z 0 N Y H A 4 / 3 Z p i Z F y S C a + O 6 3 0 5 h b X 1 j c 6 u 4 X d r Z 3 d s / K B 8 e t X S c K o Z N F o t Y d Q K q U X C J T c O N w E 6 i k E a B w H Y w v p v 5 7 S d U m s f y w U w S 9 C M 6 l D z k j B o r N W 7 6 5 Y p b d e c g q 8 T L S Q V y 1 P v l r 9 4 g Z m m E 0 j B B t e 5 6 b m L 8 j C r D m c B p q Z d q T C g b 0 y F 2 L Z U 0 Q u 1 n 8 0 O n 5 M w q A x L G y p Y 0 Z K 7 + n s h o p P U k C m x n R M 1 I L 3 s z 8 T + v m 5 r w 2 s + 4 T F K D k i 0 W h a k g J i a z r 8 m A K 2 R G T C y h T H F 7 K 2 E j q i g z N p u S D c F b f n m V t C 6 q n l v 1 G p e V 2 m 0 e R x F O 4 B T O w Y M r q M E 9 1 K E J D B C e 4 R X e n E f n x X l 3 P h a t B S e f O Y Y / c D 5 / A I z N j M E = < / l a t e x i t > = < l a t e x i t s h a 1 _ b a s e 6 4 = " D m N M k 2 Y X Y y T v + z N 0 r W X Y w O / Y u L g = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 I t Q 9 O K x B f s B b S i b 7 a R d u 9 m E 3 Y 1 Q Q n + B F w + K e P U n e f P f u G 1 z 0 N Y H A 4 / 3 Z p i Z F y S C a + O 6 3 0 5 h b X 1 j c 6 u 4 X d r Z 3 d s / K B 8 e t X S c K o Z N F o t Y d Q K q U X C J T c O N w E 6 i k E a B w H Y w v p v 5 7 S d U m s f y w U w S 9 C M 6 l D z k j B o r N W 7 6 5 Y p b d e c g q 8 T L S Q V y 1 P v l r 9 4 g Z m m E 0 j B B t e 5 6 b m L 8 j C r D m c B p q Z d q T C g b 0 y F 2 L Z U 0 Q u 1 n 8 0 O n 5 M w q A x L G y p Y 0 Z K 7 + n s h o p P U k C m x n R M 1 I L 3 s z 8 T + v m 5 r w 2 s + 4 T F K D k i 0 W h a k g J i a z r 8 m A K 2 R G T C y h T H F 7 K 2 E j q i g z N p u S D c F b f n m V t C 6 q n l v 1 G p e V 2 m 0 e R x F O 4 B T O w Y M r q M E 9 1 K E J D B C e 4 R X e n E f n x X l 3 P h a t B S e f O Y Y / c D 5 / A I z N j M E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D m N M k 2 Y X Y y T v + z N 0 r W X Y w O / Y u L g = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 I t Q 9 O K x B f s B b S i b 7 a R d u 9 m E 3 Y 1 Q Q n + B F w + K e P U n e f P f u G 1 z 0 N Y H A 4 / 3 Z p i Z F y S C a + O 6 3 0 5 h b X 1 j c 6 u 4 X d r Z 3 d s / K B 8 e t X S c K o Z N F o t Y d Q K q U X C J T c O N w E 6 i k E a B w H Y w v p v 5 7 S d U m s f y w U w S 9 C M 6 l D z k j B o r N W 7 6 5 Y p b d e c g q 8 T L S Q V y 1 P v l r 9 4 g Z m m E 0 j B B t e 5 6 b m L 8 j C r D m c B p q Z d q T C g b 0 y F 2 L Z U 0 Q u 1 n 8 0 O n 5 M w q A x L G y p Y 0 Z K 7 + n s h o p P U k C m x n R M 1 I L 3 s z 8 T + v m 5 r w 2 s + 4 T F K D k i 0 W h a k g J i a z r 8 m A K 2 R G T C y h T H F 7 K 2 E j q i g z N p u S D c F b f n m V t C 6 q n l v 1 G p e V 2 m 0 e R x F O 4 B T O w Y M r q M E 9 1 K E J D B C e 4 R X e n E f n x X l 3 P h a t B S e f O Y Y / c D 5 / A I z N j M E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D m N M k 2 Y X Y y T v + z N 0 r W X Y w O / Y u L g = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 I t Q 9 O K x B f s B b S i b 7 a R d u 9 m E 3 Y 1 Q Q n + B F w + K e P U n e f P f u G 1 z 0 N Y H A 4 / 3 Z p i Z F y S C a + O 6 3 0 5 h b X 1 j c 6 u 4 X d r Z 3 d s / K B 8 e t X S c K o Z N F o t Y d Q K q U X C J T c O N w E 6 i k E a B w H Y w v p v 5 7 S d U m s f y w U w S 9 C M 6 l D z k j B o r N W 7 6 5 Y p b d e c g q 8 T L S Q V y 1 P v l r 9 4 g Z m m E 0 j B B t e 5 6 b m L 8 j C r D m c B p q Z d q T C g b 0 y F 2 L Z U 0 Q u 1 n 8 0 O n 5 M w q A x L G y p Y 0 Z K 7 + n s h o p P U k C m x n R M 1 I L 3 s z 8 T + v m 5 r w 2 s + 4 T F K D k i 0 W h a k g J i a z r 8 m A K 2 R G T C y h T H F 7 K 2 E j q i g z N p u S D c F b f n m V t C 6 q n l v 1 G p e V 2 m 0 e R x F O 4 B T O w Y M r q M E 9 1 K E J D B C e 4 R X e n E f n x X l 3 P h a t B S e f O Y Y / c D 5 / A I z N j M E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D m N M k 2 Y X Y y T v + z N 0 r W X Y w O / Y u L g = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 I t Q 9 O K x B f s B b S i b 7 a R d u 9 m E 3 Y 1 Q Q n + B F w + K e P U n e f P f u G 1 z 0 N Y H A 4 / 3 Z p i Z F y S C a + O 6 3 0 5 h b X 1 j c 6 u 4 X d r Z 3 d s / K B 8 e t X S c K o Z N F o t Y d Q K q U X C J T c O N w E 6 i k E a B w H Y w v p v 5 7 S d U m s f y w U w S 9 C M 6 l D z k j B o r N W 7 6 5 Y p b d e c g q 8 T L S Q V y 1 P v l r 9 4 g Z m m E 0 j B B t e 5 6 b m L 8 j C r D m c B p q Z d q T C g b 0 y F 2 L Z U 0 Q u 1 n 8 0 O n 5 M w q A x L G y p Y 0 Z K 7 + n s h o p P U k C m x n R M 1 I L 3 s z 8 T + v m 5 r w 2 s + 4 T F K D k i 0 W h a k g J i a z r 8 m A K 2 R G T C y h T H F 7 K 2 E j q i g z N p u S D c F b f n m V t C 6 q n l v 1 G p e V 2 m 0 e R x F O 4 B T O w Y M r q M E 9 1 K E J D B C e 4 R X e n E f n x X l 3 P h a t B S e f O Y Y / c D 5 / A I z N j M E = < / l a t e x i t > + + · · · + + · · · + iM 3 = iD + iM 3,df < l a t e x i t s h a 1 _ b a s e 6 4 = " l u T h B b X C J Z 7 H 0 I X Y S 4 K N J X D L t 1 E = " > A A A U l 3 i c v Z h b b 9 s 2 F M f d 7 t Z m t 3 T N y z B g E J Y E 6 D r X s J N s 3 Q Y M a 5 y r 2 9 y d W x M F A S V R F h d d H J J K n A p 6 2 6 f Z 6 / Z l 9 m 1 2 K D s j R c p B 9 j I B i W X + / u Q h D 8 l D H j v 9 k D D e b P 7 9 4 O F 7 7 3 / w 4 U e P H k 9 8 / M m n n 3 0 + + e S L Q 5 a k 1 M U H b h I m 9 N h B D I c k x g e c 8 B A f 9 y l G k R P i I + d i S f C j K 0 w Z S e J 9 f t P H Z x H q x c Q n L u J Q d D 7 5 N b E 3 3 f N s P r d + s Y i 9 7 F r f W a O S u u 3 5 + f n k d L P R L B 7 L f G m N X q Z r o 2 f n / M l T Z H u J m 0 Y 4 5 m 6 I G D t t N f v 8 L E O U E z f E + Y S d M t x H 7 g X q 4 V N 4 j V G E 2 V l W D C S 3 Z q H E s / y E w l / M r a J U r Z G h i L G b y A F l h H j A d C Y K q 9 h p y v 0 f z z I S 9 1 O O Y 3 d o y E 9 D i y e W 8 I r l E Y p d H t 7 A C 3 I p g b 5 a b o A o c j n 4 b q J k J g A 9 p d i v s 9 T x S S + l u D 5 g A M t j y y K X c N w P U 1 Y v R l E X f e J J E r J 6 j 6 J + Q N x B W e 9 Q d I F 5 3 X G i O g O / B d j T 2 k t D T m h y X f Y g j n u w S I K z z E E O D r U a Y J E y n 5 V L c Q x z Q x H X u l s 1 h M p + X E P v Y d h a N w b 9 h I p Z 9 n 5 L G X c S G B s 8 s / d 4 J m Y f d w O o 6 6 a c 3 a / C x K z V 4 S g k L s x w j K / d J I p Q 7 G U 2 a e c Z / O d W O 8 8 1 s j w i y w Z 5 M y J v D N I d k a 5 B 9 k d k H w h 0 p p 2 E n g + + s y Z s i l U h y j N b z I H j Z 0 g 0 U q a O p I 5 J X U l d k 3 q S e i a 9 k P R C 7 3 x f s r 7 O L i W 7 N F u l k l K 9 5 k C y g c 5 u J L v R 2 T v J 3 u l s U b J F n b U l M y Z 7 S b I l n b 2 R 7 I 0 5 w h 1 J d / S a u 5 L t 6 u x Y s m O d v Z X s r c 5 O J D v R W V O y p s 5 a k r V 0 N i f Z n M 7 m J Z v X 2 Y J k C z q L J Y t 1 F k k W j b Z D C D H B S R D 1 i p 1 R l i c x v q 3 g R G b v F x 1 H w Y v m F L V L A j H 3 s + X J L 3 F j / v d K e M 9 Y O n f X 3 i r h L W O J l L C x S k 5 K + G T o r C W I Y 6 P j w N J c c R s B X B R W b A C V m l t A p c Y w l l V q R M Q V l a 7 o d F W l q z p d U + m a T t d V u q 7 T j k o 7 O n 2 t 0 t f G t l a p E c k 3 V L q h 0 0 2 V b h p T r l J j x r d V u q 3 T H Z W a A U W l x m L Z U 6 m x T r s q N c 6 n f Z X u 6 / R A p Q c 6 P V T p o U 6 P V H p k h E G V m o F Q p W Y o V G m x N 2 B z b M J 3 a 7 s v L i s J X M L s Z Q w X S o p F 8 W 0 p H E s Q U P b E v a W S E o h P n W g c 5 X C m 7 d N q + h z c n E a d m O f Z 7 I R l W c P e B Q l x s f g O o S w R e z f O Z m y P s H 6 I b h i / C b H N 0 m j G d q k d E A 9 f E 4 8 H G i c x n 5 E Q l D B Y t T m b I s I w 3 K C y R m v B D j D p B d A F m 4 F j h q U v c z D J 8 Y B L e 3 m u W Z T 4 P 5 i b q 7 L 2 g 7 D G 4 E b c 1 + 2 V i w s 7 / 0 P r 2 r p B D s t P W 2 e Z j W M G F 3 E x R 5 k d Q r 7 D p 1 s 2 F Z / 6 S r u t 4 M D Z J D K J J M y s 6 Z Z l n G u O P 9 S N D r c q T X I 1 l C R g R 6 R f l a J r P h S B / z E n o T d G F U h V g H i l x h l p o D e 2 1 G g i R i J 8 + Q 7 T x H T L g A r / I w o 5 x O m 0 g E F x 5 X + e N b 5 3 4 e y 2 m r n 9 b x J Q a p K N X D a a E T k v t 9 0 U m 3 W N J m m f W V D D W g x 7 G L I Z 7 X r O U n 6 d w H z D m 2 9 1 D / J n c 9 / O 6 J Z S S F g w V k X z F a I t V b B l C N C t p W L u f E i r M p Z W m U P S X l l q G k V D q 2 V Z Y R q G v g G p d G w h x 6 H 4 i h R p t T Z 0 F z I B F z W q 3 e v 6 A P 0 x E P f y D D d 6 4 y i H s A n / x m A C I y M N P L Y y C k V t C 4 V S A a N Z G U A E J C J 9 1 0 Y B G S a C g N l G b U S r m 4 Q M F D L 8 r C 0 + x i l Y p w O K l W 6 n c 4 f i V j J G I 3 5 X y L O l j c V u a 2 6 M I o n 6 i A n R 9 u b O Y r d b r e q F K Y b U Z Q 0 + j q s V Y Q A e 3 l h f G k u d A j u q A z c J c 8 s 6 D 6 a 4 i M k 0 y u B d v 1 6 L + 3 W c Q F r u Y K o x j / h q X f N C z A O R j I 4 4 f D H i w T U s 2 + 7 M C z 0 N F M U 7 R r E n i p e N Y l 8 U r x r F D h Y Z o 4 N 7 J M 7 w Z V q s f N 0 + L j Q Y X s c p 4 o R j V n i A x B 6 s O h H g r K 2 i 8 G f w Z 0 l 7 g a / y U l A b x m c Y u l 2 / w I B p b L 1 o N l o 4 s u D y M i a e R X c 0 s n n f R n p 3 N L J 2 3 0 b 4 H Y 3 s 3 9 H I + e R 0 S / 8 J z n w 5 n I M b R O O n 3 Y X p V + 3 R z 3 O P a l / V v q k 9 q 7 V q L 2 u v a u u 1 n d p B z a 3 9 X v u j 9 m f t r 6 k v p 3 6 d W p 1 a H 0 o f P h j V e V o r P V O 7 / w D 7 M c e 5 < / l a t e x i t > iM 3 = iD + iM 3,df < l a t e x i t s h a 1 _ b a s e 6 4 = " l u T h B b X C J Z 7 H 0 I X Y S 4 K N J X D L t 1 E = " > A A A U l 3 i c v Z h b b 9 s 2 F M f d 7 t Z m t 3 T N y z B g E J Y E 6 D r X s J N s 3 Q Y M a 5 y r 2 9 y d W x M F A S V R F h d d H J J K n A p 6 2 6 f Z 6 / Z l 9 m 1 2 K D s j R c p B 9 j I B i W X + / u Q h D 8 l D H j v 9 k D D e b P 7 9 4 O F 7 7 3 / w 4 U e P H k 9 8 / M m n n 3 0 + + e S L Q 5 a k 1 M U H b h I m 9 N h B D I c k x g e c 8 B A f 9 y l G k R P i I + d i S f C j K 0 w Z S e J 9 f t P H Z x H q x c Q n L u J Q d D 7 5 N b E 3 3 f N s P r d + s Y i 9 7 F r f W a O S u u 3 5 + f n k d L P R L B 7 L f G m N X q Z r o 2 f n / M l T Z H u J m 0 Y 4 5 m 6 I G D t t N f v 8 L E O U E z f E + Y S d M t x H 7 g X q 4 V N 4 j V G E 2 V l W D C S 3 Z q H E s / y E w l / M r a J U r Z G h i L G b y A F l h H j A d C Y K q 9 h p y v 0 f z z I S 9 1 O O Y 3 d o y E 9 D i y e W 8 I r l E Y p d H t 7 A C 3 I p g b 5 a b o A o c j n 4 b q J k J g A 9 p d i v s 9 T x S S + l u D 5 g A M t j y y K X c N w P U 1 Y v R l E X f e J J E r J 6 j 6 J + Q N x B W e 9 Q d I F 5 3 X G i O g O / B d j T 2 k t D T m h y X f Y g j n u w S I K z z E E O D r U a Y J E y n 5 V L c Q x z Q x H X u l s 1 h M p + X E P v Y d h a N w b 9 h I p Z 9 n 5 L G X c S G B s 8 s / d 4 J m Y f d w O o 6 6 a c 3 a / C x K z V 4 S g k L s x w j K / d J I p Q 7 G U 2 a e c Z / O d W O 8 8 1 s j w i y w Z 5 M y J v D N I d k a 5 B 9 k d k H w h 0 p p 2 E n g + + s y Z s i l U h y j N b z I H j Z 0 g 0 U q a O p I 5 J X U l d k 3 q S e i a 9 k P R C 7 3 x f s r 7 O L i W 7 N F u l k l K 9 5 k C y g c 5 u J L v R 2 T v J 3 u l s U b J F n b U l M y Z 7 S b I l n b 2 R 7 I 0 5 w h 1 J d / S a u 5 L t 6 u x Y s m O d v Z X s r c 5 O J D v R W V O y p s 5 a k r V 0 N i f Z n M 7 m J Z v X 2 Y J k C z q L J Y t 1 F k k W j b Z D C D H B S R D 1 i p 1 R l i c x v q 3 g R G b v F x 1 H w Y v m F L V L A j H 3 s + X J L 3 F j / v d K e M 9 Y O n f X 3 i r h L W O J l L C x S k 5 K + G T o r C W I Y 6 P j w N J c c R s B X B R W b A C V m l t A p c Y w l l V q R M Q V l a 7 o d F W l q z p d U + m a T t d V u q 7 T j k o 7 O n 2 t 0 t f G t l a p E c k 3 V L q h 0 0 2 V b h p T r l J j x r d V u q 3 T H Z W a A U W l x m L Z U 6 m x T r s q N c 6 n f Z X u 6 / R A p Q c 6 P V T p o U 6 P V H p k h E G V m o F Q p W Y o V G m x N 2 B z b M J 3 a 7 s v L i s J X M L s Z Q w X S o p F 8 W 0 p H E s Q U P b E v a W S E o h P n W g c 5 X C m 7 d N q + h z c n E a d m O f Z 7 I R l W c P e B Q l x s f g O o S w R e z f O Z m y P s H 6 I b h i / C b H N 0 m j G d q k d E A 9 f E 4 8 H G i c x n 5 E Q l D B Y t T m b I s I w 3 K C y R m v B D j D p B d A F m 4 F j h q U v c z D J 8 Y B L e 3 m u W Z T 4 P 5 i b q 7 L 2 g 7 D G 4 E b c 1 + 2 V i w s 7 / 0 P r 2 r p B D s t P W 2 e Z j W M G F 3 E x R 5 k d Q r 7 D p 1 s 2 F Z / 6 S r u t 4 M D Z J D K J J M y s 6 Z Z l n G u O P 9 S N D r c q T X I 1 l C R g R 6 R f l a J r P h S B / z E n o T d G F U h V g H i l x h l p o D e 2 1 G g i R i J 8 + Q 7 T x H T L g A r / I w o 5 x O m 0 g E F x 5 X + e N b 5 3 4 e y 2 m r n 9 b x J Q a p K N X D a a E T k v t 9 0 U m 3 W N J m m f W V D D W g x 7 G L I Z 7 X r O U n 6 d w H z D m 2 9 1 D / J n c 9 via a sequence of one-particle exchanges and M (u,u) df,3 is driven by the three-body K matrix K df,3 representing shortdistance three-particle interactions. Here k and p are the momenta of one of the hadrons in the initial and final state, respectively. We refer to this hadron as the spectator. The other two hadrons, called a pair, associated with the given spectator are projected to definite angular momentum. Since only total angular momentum is conserved, the unsymmetrized amplitude is a non-diagonal matrix in the pair's angular momentum space. In addition to the external momentum and angular dependencies, the amplitude depends on the total center-of-momentum (CM) energy of the three-particle system, denoted by E, which is suppressed in the argument list of Eq. (1).
From here on we make two simplifications. First, we assume that the two-body subsystem contains contributions from the = 0 partial wave only, leaving just one matrix element of Eq. (1) for our consideration. The method of the numerical solution presented here applies to any partial wave amplitude, but we constrain ourselves to the S wave case just for simplicity. Although here we are primarily interested in the unsymmetrized amplitude, the fully symmetrized amplitude is obtained by summing over the nine possible spectator momenta, which for pairs in S wave reduces to where P p = {p, a , −p − a } and P k = {k, a, −k − a} and a and a are the momenta of one of the particle in the initial and final pair state, respectively. The second simplification we describe below. The first term appearing in Eq. (1) is D, which represents the sum over all possible pair-wise interactions mediated by one-particle exchanges, is defined by the integral equation (see Fig. 1), Here M 2 is the S wave 2 → 2 scattering amplitude describing the initial and final interactions among particles in the pair. Their energy is fixed by the momentum of the spectator, where ω k ≡ √ m 2 + k 2 , k ≡ |k|, and the E 2,k is the energy of the pair evaluated in their CM frame. The ladder amplitude D is driven by the exchange propagator G, which describes the long-range interactions between the intermediate pair and spectator, and is defined as and H(p, k) is a cut-off function to render the integral in Eq. (3) finite. We use the cut-off defined in Ref. [46], such that the function evaluates to unity in the physical region and smoothly interpolates to zero at E 2,k = 0. It is advantageous to work with the amputated amplitude d, which reduces the singularities associated with the initial and final scattering of the two-particle system associated with the pair. Following Ref. [56], we define it as: It is natural to expect that it is less susceptible to instabilities when evaluated numerically close to poles of the twobody amplitudes M 2 . It is important to emphasize that this is not an approximation, but a definition of d, which is more advantageous for systems where the two-body subsystems have either bound states or resonances. Using above equation and Eq. (3), we see that d satisfies Equation (7) makes evident the claim that d is less sensitive to singularities associated with the initial/final twoparticle states. In particular, if these couple to bound states, D will have poles on the real axis while d will not. Nevertheless, d does depend on M 2 , and as a result d can still exhibit singular behavior at the two-particle thresholds. In Sec. VI we provide numerical evidence of the manifestation of these singularities. Furthermore, we give illustrative comparisons between d and D for kinematics near the two-particle thresholds and a two-body bound state.
The ladder amplitude does not contain any information about short-range three-body physics. This is described by the amplitude M df,3 , which is the second term of Eq. (1). The short-distance interactions are encoded into a three body K-matrix, denoted K df,3 , which is the driving term for the integral equation for M df,3 . This equation depends on K df,3 as well as the ladder amplitude D (see for example the discussion in Sec. V of Ref. [46]). Therefore, within the framework of Ref. [46] one must determine D first, and then for a given K df,3 the second integral equation for M df,3 can be solved.
This bring us to the second simplification we make in this study. From this point forward, we assume that the threebody K matrix is zero, thus the scattering amplitude is dominated by exchanges between two-particle subprocesses. Given this assumption, the explicit dependence of M df,3 on K df,3 is not needed here, only the fact that as K df,3 → 0, also M df,3 → 0, and therefore the three body amplitude is reduced to its ladder part: Our focus here is to develop an efficient framework for evaluating the integral equation for the ladder amplitude in the most singular scenario, namely when the two-body amplitude M 2 has a bound state pole present in the integration range of Eq. (3). As we have already emphasized, it is more efficient to determine the amputated amplitude of Eq. (7), as in this case the external pole contributions are absent from the calculation, and thus in the remainder of this work we focus on d. Having determined d, including a non-zero K df,3 contribution is straightforward. This would, of course, require first determining K df,3 either from the lattice QCD spectrum [45] or experimental data. To summarize, we consider two key approximations. The two-body subsystem is saturated by = 0 and the K df,3 = 0. This is the same limit considered by Ref. [75]. As a result, in Sec. III, we are able to provide direct comparison of the numerical solutions of the integral equations obtained here with those obtained in the aforementioned reference. The advantage of the techniques presented here are multifold. First, the numerical solutions can be reached with timeframes that are orders of magnitude shorter. This allows to map these functions continuously. Second, the solutions are systematically improvable. Third, the framework presented here holds for any kinematics above the three-particle thresholds, which is not the case for Ref. [75].
A. The J = 0 scattering amplitude The integral equation for D (or for equivalently d) is taxing due to two main issues. First, the integrand is singular. This can be seen, for example, in Eq. (8) where the two-body amplitude M 2 appears under the integral, which in general can have branch cuts and poles. Secondly, the integral is three-dimensional, which makes the problem of the numerical solution of the equation much more complex.
We can avoid the second complication by employing the partial wave projection in total angular momentum J, which leads to an infinite number of one-dimensional integral equations that are simpler to evaluate. In practice for a system with definite quantum numbers, only a finite number of these equations must be solved. Since here we consider the limit where the two-particle subsystem has a single non-zero partial wave, namely = 0, the only source of angular dependence arises from the relative momentum between the two-particle subsystem and the spectator. In the CM frame, it is given by an angle of the spectator momentum. In other words, where m J is the projection onto some external z axis. From angular momentum conservation, the resultant amplitude must be diagonal in J and from azimuthal symmetry it should be independent of m J , For simplicity, we only consider the J = 0 component of the total amplitude, which is denoted by a subscript S (for S-wave). With this, we can define the J = 0 component of D by integrating out the overall angle dependence of Eq. (3) where we have introduced We used the fact that H is independent of the angle, see Eq. (5). By partial-wave projecting the exchange propagator, we have effectively softened its singularity from a pole to a logarithm. As one would expect, this further simplifies the numerical evaluation of the integral equations. We remark that having a non-zero value of , was necessary to define the integral in Eq. (13). As we will see below, having a non-zero value of will play an important role in obtaining numerical solutions of the integral equations. Ultimately the desired solutions can be obtained by taking the → 0 limit of the subsequent solutions.
Since the initial and final state two-particle scattering amplitudes are unaffected by the angular projection performed above, the function d defined by Eq. (7) has an expansion similar to Eq. (10), and the integral equation for the partial wave projected d, Eq. (8), reads In the remainder of this work, we consider this form of the ladder equation. Equation Eq. (12) was solved in Ref. [69] to obtain the first three-body amplitudes from lattice QCD. Formally, successive iterations of Eq. (3) yield a solution in terms of an increasing number of exchanges between the two-particle subsystems. For weakly coupled two-body systems, i.e. for small ma, the series rapidly converges and the first few orders dominate the solution. In the case considered here, i.e. for strongly interacting systems in which the twoparticle subsystem forms a bound state, the perturbation series fails to converge and we are forced to resort to a non-perturbative numerical approach.

B. Integral equations in the presence of a two-body bound state
As announced in Sec. I, here we are interested in the implications of these integral equations for systems where the two-body subsystem can become bound. In this work, we consider the effective range expansion for M 2 , and assume it is dominated by the leading order (LO) order term where q 2k = E 2 2,k /4 − m 2 is the relative momentum between the two particles in their CM frame and a is the scattering length. Figure 2 shows plots of |M 2 | as a function of E 2k /m = √ s 2k /m for ma = 2, 6, and 16, which are the three cases we study in detail in the subsequent solutions of the integral equations. The two-body scattering amplitude has a pole on the real s 2k axis, which we call s b . Near the bound state, the scattering amplitude takes the form, where g is the residue at the pole, which can be interpreted as the bound state wave function renormalization factor. Since we are considering the contribution of a pole in the first Riemann sheet below the threshold, the relative momentum of the two-particle subsystem lies on the positive imaginary axis, q 2k = iκ 2k , where κ 2k > 0 is the binding momentum. In terms of the LO effective range expansion Eq. (15), the binding momentum is κ 2k = 1/a. In general, the pole position in s 2k can be written in terms of the binding momentum in the standard way, By equating Eq. (15) to Eq. (16), one finds that the residue in the LO in the effective range expansion is given by: The integral equations of interest are written in terms of the spectator momentum, therefore we need to define the value of k corresponding to the bound state pole. We will label this "on-shell" value of k as q. This can be obtained by fixing the two-particle subsystem to be at the bound state pole, and by requiring that systems to be in its total CM frame, the sum of the energy of the two-body subsystem and the spectator satisfies E = s b + q 2 + m 2 + q 2 . Solving for q gives: where λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + yz + zx) is the Källén triangle function. Consequently, the integral equation presented in Eq. (14) involves an integral over this pole, which makes the numerical convergence of the solutions harder to achieve. As discussed above in the context of partial wave projection [see Eq. (13)], it was necessary to introduce a non-zero value of . Here we once again are required to introduce a non-zero value of to mitigate the real-axis pole of M 2 . This shifts the pole slightly away from the axis of integration, allowing for a rigorous definition of the integral. Below the three-particle threshold the exchange propagator never goes on-shell, and thus is a smooth function for all energies in this domain, leaving the two-particle bound state pole as the only singularity in the integration region. Given that it is one of the main issues encountered and addressed in this work, from here on we make the dependence explicit in the quantities that are most sensitive to its presence. To that end, we shall denote the S wave projected exchange propagator by G S (p, k) → G S (p, k; ), where G S (p, k; ) is defined exactly as in Eq. (13). The shift is also implemented in the energy of the two-body system in the following way: This shift propagates through to the amplitudes D and d, where the latter's integral equation is given by where d (u,u) S (p, k) is given by the limit as → 0. After solutions of the integral equation are obtained, it is necessary to analyze their → 0 limit numerically.
Given the definition of D S in Eq. (12) one can see that the three-body scattering amplitude has external poles associated with the two-body states. The residue of these is related to the scattering amplitude between the spectator and the two-body bound state. Labeling the spectator as "ϕ" and the bound state as "b", we label this amplitude as M ϕb , and denote it as the 2+1 scattering amplitude. More specifically, by continuing the initial and final two-particle subsystems to the bound state poles, the three-body scattering amplitude is related to M ϕb by the LSZ reduction, This implies that by evaluating the integral equation for the M (u,u) 3,S as a function of energy one will see this double-pole structure. By determining the residues of these, using the numerical equivalent of one can determine the ϕb → ϕb scattering amplitude. One important point is the fact that these identities concern the unsymmetrized amplitude. One can obtain the symmetrized amplitude following the procedure defined in Eq. (2), with eight of the nine terms not contributing to the final sum. This procedure is well defined, however, it requires an additional step in the numerical determination of the solutions of the singular integral equations. Namely, one needs to scan M (u,u) 3,S as a function of the two-body energy and obtain the residue at the bound state pole. This might be susceptible to numerical instabilities since at this point the zeroes in Eq. (23) need to cancel the divergent pole terms. Fortunately, this numerical issue disappears when performing a redefinition of the three-particle amplitude, as in Eq. (7) for D. Such a redefinition allows one to evaluate analytically the cancellation of the external two-body poles and zeros. Given a numerical solution of the integral equations, one would still need to numerically evaluate the limit to the bound state pole, but this would instead be done for a smooth function. Using the expressions above and Eq. (7), it is easy to then see that in this limit, M ϕb can be obtained from d via The ϕb → ϕb amplitude is our primary subject of study, which, according to the above equation, for chosen a is completely determined by the solution of Eq. (14). By construction, solutions of Eq. (14) satisfy the three-particle S matrix unitarity [51]. Below the three-particle threshold, the ϕb → ϕb amplitude in turn satisfies the standard 2 → 2 S matrix unitarity, where is the two-body phase space between the bound state and the spectator. It follows from Eq. (25) that the amplitude is bounded by unity as |ρ ϕb (E)M ϕb (E)| ≤ 1 in this kinematic region. Additionally, Eq. (25) tells us that the amplitude can be written using a 2 → 2 K-matrix construction, where Re M −1 ϕb (E) ≡ K −1 ϕb (E) is real below the 3ϕ threshold. The ϕb → ϕb K matrix can be written in terms of a real phase shift δ ϕb in the standard way At the ϕb threshold, q cot δ ϕb reduces to the bound state-spectator scattering length, which we label as b 0 , In the Sec. III, we present numerical results for the ϕb → ϕb amplitude below the three-particle threshold. We compare our findings to Ref. [75], which studies the same ϕb scattering system from the perspective of the associated finite volume formalism [45].

III. NUMERICAL RESULTS
Having identified the final form of the integral equation we wish to evaluate, namely Eq. (14) for the partial-wave projected, amputated amplitude d S , we proceed to take the same first steps as in Ref. [69]. The first step amounts to discretize the momenta appearing in the integral equation. This allows one to write the integral equation as a matrix equation that could be solved numerically. We convert the integral to a sum over equidistant mesh points, where the number of points being denoted by N . For each value of E/m, we compute the solution for various values of N , ranging from 1000 to 6000. The solutions of the original integral equations are recovered by extrapolating N → ∞ while keeping N fixed. This assures that results are insensitive to both N and . As one can expect, the results may converge faster or slower depending on the points chosen in the (N, ) plane. In Sec. V, we explain the trajectories chosen and the reasoning behind them. We study three cases for the two-particle scattering amplitude, namely ma = 2, 6, and 16. We first summarize how solutions are computed, and relegate details for the interested reader to the following sections.
Because of the presence of the pole in M 2 , to gain confidence we carry out this procedure in two different ways. In the first method, we follow the procedure outlined in the previous section, where the bound state pole is moved off the   real axis. We refer to this as the "brute force" (BF) method. In the second procedure, we evaluate the contribution of the pole analytically. This results in a modified integral equation with a less singular kernel, that we then solve numerically. We refer to this as the "semi-analytic" (SA) method. Detailed descriptions of both of these methods are given in Sec. IV. As discussed in the previous section, solutions must satisfy the S matrix unitarity Eq. (25). We use this fact as a check on the quality of solutions as a function of N . For each value of E, the deviation away from this condition gives us an estimate of the systematic error. In particular, we define which gives a percent measure of the deviation of the solution from unitarity. As ∆ρ ϕb → 0, the solution better satisfies the unitarity relation. Therefore, we use this measure to scan for satisfactory solutions and improve the candidate solution by varying and N parameters to drive ∆ρ ϕb as small as possible. Our goal in this study is to achieve sub-percent level deviation, on the order of O(10 −1 ) − O(10 −3 )%. For most energies, this is easily attainable, with exceptions occurring around points where either ρ ϕb or the amplitude vanish.
Since the result of the matrix equation depends on N and , we must take the ordered double limit as first N → ∞ and then → 0. The order of these limits can not be reversed. Qualitatively, as steadily decreases the bound state pole moves closer to the axis of integration. This results in the higher and narrower singularity of the integrand in Eq. (14) and requires greater mesh sizes N to probe the integration kernel effectively. We express as a function of N in the form ∝ 1/N , which we derive in Sec. V A. The proportionality constant is a function of E and includes a controllable parameter which we call η. We show numerical evidence for this behavior of and show that for some restricted values of η, the resulting solutions are insensitive to our desired working precision, which is a sub-percent level deviation of unitarity. Finally, we perform large N extrapolations to estimate the N → ∞ limit. We find that these extrapolations improve the deviation from unitarity by a few orders of magnitude.
We proceed to the presentation of solutions to the integral equations for total CM energy in the domain 1+ √ s b /m < E/m < 3. The first case examined is ma = 2, which describes a deeply bound state in the two-body subsystem. The resulting bound state plus spectator scattering amplitude is shown in the top panel of Fig. 3. Shown with the vertical dashed line is the ϕb threshold at E/m = 1 + √ 3. Since below the three-particle threshold, the system is a two-particle system, it is constrained by the usual S matrix principles such as |ρ ϕb M ϕb | being bounded by unity. The corresponding q cot δ ϕb is shown in the middle panel of Fig. 3, along with points computed of the same system but using the finite volume formalism developed for three-particle scattering processes, Ref. [75]. In that work, the authors analytically continue the quantization condition to the ϕb system. Fitting our solution to an effective range expansion, where b 0 is the S-wave scattering length and r 0 is the effective range of the bound state plus spectator system, we find the fit values mb 0 ≈ 6.4 and mr 0 ≈ 2.3. For these kinematics, we find an excellent agreement with that study. The bottom panel of Fig. 3 shows the unitarity deviation, which is a measure of the systematic error arising from deviations from the unitarity condition Im M −1 ϕb = −ρ ϕb induced by considering a finite N and . This ∆ρ ϕb function shows sub-percent deviations from this condition and gives a measure of the quality of the solutions for a given energy E. For exact solutions, ∆ρ ϕb = 0. In our calculations for ma = 2, we find for all energy points except at threshold show sub-percent deviations, indicating that our numerical approximation is satisfactory for the precision we desire in this study. At the threshold, the deviation measures on the order of a percent, which is still exceedingly good for our working precision. Generally, near the threshold we find that solutions show a larger unitarity deviation than other energy points for all a, owing to the small kinematic phase space.
In the ma = 6 case, the two-body bound state moves toward the threshold, producing a shallow bound state in the two-body subprocess, which results in the ϕb → ϕb amplitude shown in Fig. 4. Noticeably in this case a zero of the amplitude is slightly below (E/m) 2 ≈ 8.86. This zero corresponds to a pole in q cot δ ϕb at this kinematic point. Using Eq. (31) to describe this case works in a very limited region close to the threshold due to this pole, which gives a scattering length mb 0 ≈ −3.6. Deviations from unitarity lie at the sub-percent level except at the threshold and the zero of the amplitude. However, such a large deviation from unitarity is not a major concern since it results in a percent-level systematic error for an observable which is equal to zero at this point.
Finally, for ma = 16 there is a very shallow bound state in the two-body channel. In Fig. 5, the top, middle, and bottom panels show the resulting amplitude, q cot δ ϕb , and ∆ρ ϕb , respectively. We find the corresponding mb 0 ≈ 150, consistent with Ref. [75]. However, we observe a significant deviation of the q cot δ ϕb as compared to Ref. [75] for energies near the three-body threshold. These are attributed to the failure of the method used in Ref. [75] near and above the three-particle region.
Given this slight deviation, it is worthwhile to summarize the method used in Ref. [75] and explain its errors in this kinematic region. There, finite-volume energies levels for these toy theories using large volumes, mL = 20 − 70, were obtained. For these volumes, some states lie below the 3ϕ threshold and above the ϕb threshold. If these are sufficiently below the 3ϕ, these can be approximated as two-body finite-volume states and must therefore satisfy the two-body quantization condition [12,13], from which one can obtain q cot δ ϕb . As the finite-volume approaches the 3ϕ threshold, it no longer satisfies the two-body quantization condition. In general, this transition is not an abrupt, but rather a continuous behavior. This is consistent with the deviation between the results found here, which are not susceptible to these issues, and that obtained ibid.
As a final comparison to Ref. [75], we solve the amplitude at threshold energy E = m + √ s b , and compute the ratio of the ϕb scattering length to the two-particle scattering length b 0 /a as a function of ma. This result is shown in Fig. 6, where the blue line is the result of our calculation using the extrapolation method as described above, and the open black circles are the results from Ref. [75], which uses the finite volume formalism to extract the ϕb phase shift. We also compare our result to one computed from a non-relativistic effective field theory (NREFT) formalism [76], where we follow Ref. [75] by choosing a cutoff Λ = 0.75m.cWe find that for small ma, corresponding to highly relativistic systems, the solution differs from the NREFT considerably with the presence of an additional pole as compared to the NREFT result. Near the pole at ma ≈ 13, our solution begins to deviate from the NREFT solution. Two factors can attribute to such a discrepancy, one of which is the fact that the NREFT result is computed at a finite matrix size of N = 2000 whereas we perform an extrapolation on our solutions. The second is that the solution is scheme dependent through the parameter Λ, where we chose the value shown in order to compare to the result shown in Ref. [75].

IV. DETAILS OF NUMERICAL METHODS
In this section, we discuss the procedure of converting the integral equations to matrix equations and define both the BF and SA methods introduced in the previous section. To evaluate the integral appearing in Eq. (14), it is necessary to discretize the spectator momenta in either method. We denote these momenta by a discrete index, i.e. we make the replacement k → k n , and generate a uniform mesh of points. The minimum value that the momenta can take is k min = 0, while the maximum value corresponds to the point at which the cutoff function H, defined in Eq. (5), is zero. It vanishes when E 2,k = 0, which implies k 2 max = (E 2 − m 2 )/2E 2 . We replace the measure dk → ∆k = k max (E)/N which is the distance between mesh points for a given energy. We find that because the integrand is singular, it is necessary to use a large number of mesh points to finely sample the kernel. 2

A. Brute force method
Having a uniform mesh of points, one proceeds to discretize the momentum appearing in Eq. (14). In practice, to solve this equation it is necessary to keep N fixed as a finite parameter and test the convergence with N and . Making the N dependence of d explicit we have To solve this we write d as a matrix in the space defined by the set of {k n }, with matrix elements where B is a matrix defined as B nn = δ kn,k n + ∆k k 2 n (2π) 2 ω k n G S (k n , k n ; ) M 2 (k n ; ) .
Note, although above we assumed d is a matrix, it is sufficient to assume only p ∈ {k n } and leave k as a continuous variable. Since we are interested in the ϕb → ϕb amplitude, we choose k = q. For a large value of N , using Eq. (32), one easily interpolates to any continuous value of these momenta within the kinematically allowed values, including the on-shell point q for bound state plus spectator system, defined in Eq. (19).
In Sec. V we explain how one may assess systematic errors of solution to Eq. (32). In particular, to arrive at the solution to the integral equation one must take the ordered-double limit, In practice, by calculating the amplitude for several values of sufficiently large N and small enough one could perform a careful extrapolation of the numerical result and test the convergence. To make the extrapolation of the convergence tests systematic, in Sec. V A we determine the asymptotic behavior of the error for large N and small , finding where η ≡ 2πN q /k max and q is function of the energy and is linearly proportional to defined in Eq. (49). This η parameter provides a relation between and N , which for a fixed η the ordered double limit is ensured. This explains explicitly that for a given value of N one cannot make arbitrarily small, otherwise, the error introduced will no longer be exponentially suppressed. A sufficient condition is that for a given q the matrix size N must be large enough, satisfying inequality η 1. In Sec. V C we provide numerical evidence of this behavior of the error. where is an off-shell extension of the K matrix of the ϕb system. In other words, when k n , k n = q this coincides with the physical ϕb K matrix as defined in Eq. (27). In principle, K ϕb must be a real function below the three-particle threshold. However, since we work with finite N , K ϕb may be complex. The level of complexity is measured by the same ∆ρ ϕb as defined in Eq. (30). As in the case with the BF method, one must then take the ordered, double limit defined Eq. (36). In Sec. III we showed results for the ϕb amplitude as computed using the SA method. Both the BF and SA methods reliably recover the amplitude at the various N we tested which were consistent with each other. We found that the SA method for most kinematics can yield results that have a unitarity deviation an order of magnitude lower than the BF method. Since the results at the precision at which we show plots are visually indistinguishable, we show only the results for the SA method throughout this article.

V. ASSESSING SYSTEMATICS
Since the integral equations in both methods are always solved numerically at a finite N , and for a nonzero , our solutions deviate systematically from the result defined by the N → ∞ and → 0 limits. In this section, we present a detailed discussion of these systematic effects.
A. Proof of the O e −η systematic error Our first step is to evaluate the difference between the solution of the desired integral equation, given in Eq. (14), and the numerical solution, which satisfies Eq. (32), and is obtained for a finite value of N . Both the integral equation and the matrix equation require the introduction of an , which must be set to zero at the end of the computation. Here, we consider the difference of these at the stage where is nonzero. Assuming our isotropic mesh, we find that it is given by In the second line, we used d ) and kept only the leading term. This difference is then written as a combined summation and integration using the Poisson summation formula in the third line. In general, the integral is saturated by its singularities. In this case, we have the pole singularity due to M 2 , as well as logarithmic ones. Given that the former is the largest, and consequently the source of the leading error, we approximate the integral by the contribution due to the pole which is at Above the ϕb threshold, q (E) > 0, as expected. With this, we obtain the correction near the pole, which is where in the last equality we have defined η ≡ 2π q N/k max and assumed η > 1 while ignoring contributions from n > 1 that are further suppressed. This tells us the condition needed for the systematic error to be suppressed, The parameter η characterizes dependence of as a function of N . For a fixed η, the ordered double limit in and N is ensured by a single limit in N . Our derivation does not completely fix the N dependence of the function, as we have only looked at the dominating term in the series of the NLO correction. Moreover, since the N -dependent is propagated through to the amplitude in a form of an energy shift, this gives an explicit N dependence to the energy. Moreover, although the above bound is derived for the explicit pole in M 2 , and thus is applicable to the BF method, we use the same constant η trajectory in the SA method in order to ensure the double limit is properly taken.

B. Large N extrapolations
Our solutions are obtained from the matrix equations, Eqs. (32) and (43), therefore they explicitly depend on the mesh size N . Moreover, from the previous section, we assert that ∝ 1/N , which introduces additional Ndependencies in both the BF and SA methods. As discussed above, fixing η 1 assures that the leading N behavior is due to the shift in the energy rather than the discretization. This tells us that at large N , the finite N amplitude for a given E/m, ma, and η can be represented as the N → ∞ amplitude and a correction factor in 1/N , Therefore, as N → ∞ the numerical solution better approximates the continuum solution for M ϕb . Since we always work with a finite N , we employ a program which solves the integral equation for various values of N , and perform an extrapolation in N in order to estimate the N → ∞ amplitude. The extrapolation is executed by fitting the finite N amplitude for a fixed E/m, ma, and η to a function which has a polynomial dependence on 1/N . We choose two kinds of fit models: linear and quadratic in 1/N , where α and β are complex parameters for each energy. Although the extrapolated amplitude is our best estimate of the N → ∞ solution, since our calculations are always performed at finite N , there is a systematic error that propagates through to the extrapolated value. This impacts the result, which can be measured e.g. through a ∆ρ ϕb test computed with the extrapolated value of the amplitude. In Fig. 7  on the asymptotic deviation of unitarity by a few orders of magnitude as compared with the finite N values. We see evidence that the η-dependence of the solutions are mild when compared to the finite N effects.
To obtain the final amplitudes presented in Sec. III, we generate a small ensemble of extrapolated solutions by fitting both linear and quadratic models to various subsets of our finite N results. In our computations, we used 11 equidistant points in the interval 1000 ≤ N ≤ 6000. For both the linear and quadratic models, we repeat the fit by successively excluding the lowest value of N until we reach N = 4000. This results in a total of 14 fits, from which we choose the one which minimizes the extrapolated ∆ρ ϕb . To estimate a systematic error associated with the N -dependence, we take the difference between this fit and the one which maximizes ∆ρ ϕb as twice our systematic error associated with N . The impact of this error is orders of magnitude smaller than the results themselves and thus invisible in the presented results. We find, as shown in the main results, that our estimated unitarity deviation is sub-percent level for all kinematic points except for those where the amplitude is zero, and thus Eq. (30) becomes numerically ill-defined.

C. η independence of solutions
In the N → ∞ limit, any choice of η satisfying the bound Eq. (51) should converge to the same solution. In practice, we always work with finite N , and thus choices for η lead to systematic effects. To minimize these systematic deviations, we explore an example solution as a function of both and N for a given ma and E/m. In Fig. 8, we show a density plot of ∆ρ ϕb as a function of and N for both the BF and SA methods at fixed (E/m) 2 = 8 and ma = 2 for 100 ≤ N ≤ 1000, along with curves at constant η = 5, 15, and 25. For the BF method, for fixed N and small enough the ∆ρ ϕb grows suddenly, meaning that we moved past the optimal (N ) trajectory given by some optimal η, empirically verifying the approximate bound Eq. (51). Moreover, we find clear signal of oscillating solutions for η 10, which dampen as N increases.
We find that for η ∼ 10, there remains residual oscillations for the SA method as well, however the magnitude is considerably less than that of the BF and at the scale presented in Fig. 8 is indistinguishable. Figure 8 also shows that for too large η, ∆ρ ϕb increases. Since we work with systems with 1000 ≤ N ≤ 6000, we choose to compute solutions with 10 ≤ η ≤ 25 for both the BF and SA methods. We find that to our working precision, solutions computed with these η yield consistent results. Deviations for different η are an order of magnitude smaller than systematics estimated from the N dependence of the solutions, as illustrated in Fig. 7, therefore we show results with η = 15 and absorb any η fluctuations to the systematic errors arising from the N dependence. ϕb→3ϕ and d (u,u) amplitudes for fixed ma = 2, η = 15, and N = 5000 as a function of s2p/m 2 for three particle energies (E/m) 2 = 10, 15, and 20. The vertical dashed line in the top two panels shows the location of the two body bound state s b for this ma, and in all panels the two particle threshold at s2p = 4m 2 is indicated by the open circle on the real axis.

VI. ABOVE THE THREE-PARTICLE THRESHOLD
So far, we have considered energies below the three-particle breakup threshold. In this section, we show that the integral equations can be solved above the three-particle threshold, which was a limitation of the method presented in Ref. [75]. In this energy region, the bound state breakup amplitude ϕb → 3ϕ is kinematically accessible. Following the discussion in Sec. II B, one obtains the ϕb → 3ϕ amplitude via Note, the final state is composed of three particles, and as a result the unsymmetrized amplitude emerges. In order to symmetrize it, one can follow a procedure similar to that of Eq.
where only three terms are summed over, since one only needs to symmetrize the final state. Using similar arguments as in Sec. II B, we can write the ϕb → 3ϕ amplitude in terms of d As before, we solve for d (u,u) S for fixed initial k = q at some E, only now the final state momentum p is now free. As an illustration of solving the integral equations above the three particle threshold, we show in Fig. 9 the resulting M ϕb→3ϕ , the d (u,u) amplitude does not exhibit this bound state pole, as originally discussed in Sec. II. The amplitudes are plotted within the allowed integration region 0 ≤ p ≤ k max (E), which in terms of s 2p corresponds to 0 ≤ s 2p /m 2 ≤ (E/m − 1) 2 .

VII. CONCLUSIONS
In this work, we presented a numerical method for solving the three-body on-shell integral equations in the presence of the two-body bound-states. Our methodology approximates the integral equations as a system of N linear equations, which are solved by usual matrix inversion techniques. The method is systematically improvable insofar as the mesh of momentum points that are used to generate the equations can be finer sampled, leading to larger systems that better converge to the N → ∞ solution. Quantitatively, the quality of solutions is measured by computing the deviation from S matrix unitarity for each N . Empirically we find the deviation ∆ρ ϕb decreases as N increases, albeit slowly for matrix sizes of the order of 1000 − 6000. In addition to finite N effects, a convergence of solutions is affected by the presence of the two-body bound state, which produces a pole singularity in the region of integration.
We introduce two methods to circumvent the bound state pole. In the BF method, we regulate the pole singularity by introducing a finite which shifts the pole slightly off the real energy axis, and recover our solutions in the ordered limit of → 0 then N → ∞. Alternatively, in the SA method, we remove the imaginary part of the pole explicitly using an -regulated delta function. Using the NLO solution, we argued that the dominant systematic error is exponentially suppressed with η ∝ N and we empirically verified this behavior. To ensure the ordered double limit is properly taken, we computed solutions for fixed η and N . To increase the rate of convergence, and reduce the systematic error associated with finite N , we employed a strategy of extrapolating the results to the N → ∞ limit by generating a sequence of solutions with various sizes N and fitting the data to a model which is polynomial in 1/N . The extrapolated solutions drastically reduce the deviation from unitarity by several orders of magnitude as compared to the finite N results. As an assessment of the systematic error associated with this extrapolation, we generated a small ensemble of results by varying the included data in the fit. This is the dominating error of our studying, giving a sub-percent deviation from unitarity for the matrix sizes considered.
Our methodology is not limited to systems where the two-particle subsystem produces a bound state and can be adopted straightforwardly to resonating systems as well as processes with non-zero angular momenta. Additionally, the strategies presented here can be used to include the short-distance three-body K matrix, which would be determined by complementary lattice QCD calculations. In principle, one may extend these methods to complex energies to systematically search below the threshold for three-particle bound states, or on unphysical sheets for resonances. However, it is currently unknown how to consistently analytically continue the solutions outside the physical region.