Complex angular momentum, Riemann sheets, Regge poles and factorization

César Fernández-Ramírez
Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, Ciudad de México 04510, Mexico
E-mail: cesar.fernandez@nucleares.unam.mx

Prepared for the 2017 International Summer Workshop on Reaction Theory,  June 12-22, 2017, Bloomington IN, USA

Notebook definitions

Exit[];

(* Libraries *)
<<NumericalDifferentialEquationAnalysis`;
Needs["NumericalCalculus`"]

SetDirectory[NotebookDirectory[]];

(* Style options *)
StyleList = {AbsoluteThickness[2.0], AbsolutePointSize[3], #} & /@ {Blue, Red, RGBColor[0,0,0], Purple};
SetOptions[Plot, Frame -> True, Axes -> {True, False}, PlotStyle -> StyleList, AspectRatio -> 0.8, FrameStyle -> Directive[13]];

ϵ=0.000001;
epsilon=ϵ; 
(* the epsilon is needed to ensure that we are above or blow the cut when we compute the amplitudes in the physical axis *)

Kinematics and phase space

(* Threshold *)
sth[m1_,m2_]:=(m1+ m2)^2;
a[m1_,m2_]:=m1 m2/sth[m1,m2];
a0[m1_,m2_]:=Sqrt[a[m1,m2]];
 (* "Momentum" squared *)
q2[s_,m1_,m2_]:=a[m1,m2](s-sth[m1,m2]);
q2phys[s_,m1_,m2_]:=If[Re[s]>sth[m1,m2],q2[s,m1,m2],0]; 
(* "Momentum" *)
q[s_,m1_,m2_]:=Sqrt[q2[s,m1,m2]];
qphys[s_,m1_,m2_]:=Sqrt[q2phys[s,m1,m2]];
(* Angular momentum barrier *)
Bl[s_,l_,m1_,m2_]:= (q2[s,m1,m2]/(1+q2[s,m1,m2]))^l;
Blphys[s_,l_,m1_,m2_]:= (q2phys[s,m1,m2]/(1+q2phys[s,m1,m2]))^l;
(* Phase space with angular momentum barrier *)
Cl[s_,l_,m1_,m2_]:=q[s,m1,m2]Bl[s,l,m1,m2]; 
Clphys[s_,l_,m1_,m2_]:=qphys[s,m1,m2]Blphys[s,l,m1,m2];

Analytic continuation of the phase space

(* Analytically continued phase space with angular momentum barrier *)
(* L=0 *)
ρ0[s_,m1_,m2_]:=a0[m1,m2](s-sth[m1,m2])/Sqrt[-s+sth[m1,m2]];
(* L=1 *)
ρ1[s_,m1_,m2_]:=a0[m1,m2]a[m1,m2](s-sth[m1,m2])/(a0[m1,m2]+a[m1,m2]Sqrt[-s+sth[m1,m2]]);
(* L=2 *)
ρ2[s_,m1_,m2_]:=a0[m1,m2](a[m1,m2]^2)(s-sth[m1,m2])*((-s+sth[m1,m2])^1.5*(2+(1+3a[m1,m2]*(s-sth[m1,m2]))/(a[m1,m2](-s+sth[m1,m2]))^1.5))/(2(1+a[m1,m2]*(s-sth[m1,m2]))^2);
(* L=3 *)
ρ3[s_,m1_,m2_]:=a0[m1,m2](a[m1,m2]^3)(s-sth[m1,m2])((-s+sth[m1,m2])^2.5 (-8+(3(1+(10a[m1,m2](s-sth[m1,m2]))/3+5*a[m1,m2]^2 *(s-sth[m1,m2])^2))/(a[m1,m2](-s+sth[m1,m2]))^2.5))/(8(1+a[m1,m2](s-sth[m1,m2]))^3);
(* L=4 *)
ρ4[s_,m1_,m2_]:=a0[m1,m2](a[m1,m2]^4)(s-sth[m1,m2])*((-s+sth[m1,m2])^3.5(16+(5(1+(21a[m1,m2](s-sth[m1,m2]))/5+7a[m1,m2]^2 (s-sth[m1,m2])^2+7a[m1,m2]^3(s-sth[m1,m2])^3))/(a[m1,m2](-s+sth[m1,m2]))^3.5))/(16(1+a[m1,m2](s-sth[m1,m2]))^4);

(* General L phase space function *)
ρC[s_,l_,m1_,m2_]:=(1/Pi)*(a[m1,m2]^(1/2 + l)*(s - sth[m1,m2])*(Pi/((a[m1,m2] + 1/(s - sth[m1,m2]))^l*Sqrt[-s + sth[m1,m2]]) + 
    (1/((s - sth[m1,m2])*Gamma[1 + l]))*(a[m1,m2]^(-(1/2) - l)*Sqrt[Pi]*Gamma[1/2 + l]*
      ((1 + a[m1,m2]*(s - sth[m1,m2]))*Hypergeometric2F1[1, 1/2 + l, -(1/2), 1/(a[m1,m2]*(-s + sth[m1,m2]))] - 
       (3 + 2*l + a[m1,m2]*(s - sth[m1,m2]))*Hypergeometric2F1[1, 1/2 + l, 1/2, 1/(a[m1,m2]*(-s + sth[m1,m2]))]))));
(*Only valid for  (sth ≥ Re[s] || s∉Reals ) && Re[L]>-(1/2) *)

(* The one we are using *)
ρ[s_,l_,m1_,m2_]:=If[l==0||l==1||l==2||l==3||l==4,Which[l==0,ρ0[s,m1,m2],l==1,ρ1[s,m1,m2],l==2,ρ2[s,m1,m2],l==3,ρ3[s,m1,m2],l==4,ρ4[s,m1,m2]],ρC[s,l,m1,m2]];

Parameters

(* Define masses for: 1 + 2 -> 1 + 2 *)
mt1=0.5; mt2=1.0; DisplayForm[{"Threshold=", sth[mt1,mt2]}]

First, the phase space ρ(s) in the real axis

(* Define the partial wave, please try different lt *)
lt=1; 
(* Set the s limits to plot *)
 limitdown=sth[mt1,mt2];limitup =5*sth[mt1,mt2];
(* Here we plot *)
Quiet[GraphicsRow[{Plot[{Im[ρC[Res+I ϵ,lt,mt1,mt2]],Im[ρC[Res-I ϵ,lt,mt1,mt2]],Re[ρC[Res-I ϵ,lt,mt1,mt2]]},{Res,-1,limitup},PlotLabel->Style[Framed["Im ρ(s+i0)(blue),Im ρ(s-i0) (red),Re ρ(s+i0)(black)"],FontSize->14],GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}}](*See the right hand cut (blue=+i0, red=-i0) and the generated real part in black*),Plot[{Re[Cl[Res,lt,mt1,mt2]],Im[Cl[Res,lt,mt1,mt2]]},{Res,limitdown,limitup},PlotLabel->Style[Framed["Re Cl(s+i0)(blue),Im Cl(s+i0)(red)"],FontSize->14],GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}}](* Re Cl (blue) and Im Cl=0 (red) *),Plot[{Im[ρ[Res+I ϵ,lt,mt1,mt2]]-Re[Cl[Res,lt,mt1,mt2]]},{Res,limitdown,limitup},PlotLabel->Style[Framed["Im ρ(s+i0-Re Cl(s+i0)"],FontSize->14],GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}}](* Difference between the Imρ and ReCl *)},ImageSize->Full]]

(* We are computing at ±iϵ away from the real axis, that has a small influence that you see in the left figure, change ϵ and you will see the effect *)
(* The orange line marks the threshold, remember that it does not make sense to plot Cl(s) in the unphysical region *)

Try different L... Try a complex value!

lt=1+ I;
(* Set the s limits to plot *)
 limitdown=-1; limitup =5*sth[mt1,mt2];
(* Here we plot *)
Quiet[GraphicsRow[{Plot[{Re[ρC[Res-I ϵ,lt,mt1,mt2]]},{Res,limitdown,limitup},PlotLabel->Style[Framed["Re ρ(s+i0)"],FontSize->14],GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}}],Plot[{Im[ρC[Res+I ϵ,lt,mt1,mt2]]},{Res,limitdown,limitup},PlotLabel->Style[Framed["Im ρ(s+i0)"],FontSize->14],GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}}],Plot[{Abs[ρC[Res+I ϵ,lt,mt1,mt2]]},{Res,limitdown,limitup},PlotLabel->Style[Framed["Abs ρ(s+i0)"],FontSize->14],GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}}]},ImageSize->Full]]

Let’s try animation, we change L

(* Our function is analytical in both s and L complex planes*)
(* Set the s limits to plot *)
lrangedown=-0.4; (* Do not go below L=-1/2*)
lrangeup =8; (* Depends on how much computing power you have*) 
limitup =5*sth[mt1,mt2];
(* Here we plot *)
Animate[{Plot[Re[ρC[Res+I ϵ,l,mt1,mt2]],{Res,-1,limitup},PlotLabel->Style[Framed[{"Re ρ(s), l=",l}],FontSize->14],GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}},ImageSize->Medium],Plot[Im[ρC[Res+I ϵ,l,mt1,mt2]],{Res,-1,limitup},PlotLabel->Style[Framed[{"Im ρ(s), l=" ,l}],FontSize->14],GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}},ImageSize->Medium],Plot[Re[ρC[Res+I ϵ,l I,mt1,mt2]],{Res,-1,limitup},PlotLabel->Style[Framed[{"Re ρ(s), l=",l}],FontSize->14],GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}},ImageSize->Medium],Plot[Im[ρC[Res+I ϵ,l I,mt1,mt2]],{Res,-1,limitup},PlotLabel->Style[Framed[{"Im ρ(s), l=",l}],FontSize->14],GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}},ImageSize->Medium]},{l,lrangedown,lrangeup},AnimationDirection->ForwardBackward,ImageMargins->20,AnimationRate->0.5,AnimationRunning->False]

(* The vertical orange line is the threshold *)
(* First row: L is real an varies between lrangedown and lrangeup *)
(* Second: L is purely imaginary an varies between lrangedown and lrangeup *)

...then, ρ(s) in the complex s plane for fixed L

lt = 0;

(* Ranges with the plots *)
limitdownR =0; limitupR=2*sth[mt1,mt2];limitdownI=-1;limitupI=1;limitM=0;

(* We prepare the plots *)
(* Real part, 1st Riemann sheet *)
Quiet[ρfig1=Plot3D[Re[ρ[Res+ I Ims,lt,mt1,mt2]],{Res,limitdownR ,limitupR},{Ims,limitM,limitupI},ColorFunction->"Aquamarine",AxesLabel->Automatic]];
Quiet[ρfig2=Plot3D[Re[ρ[Res+ I Ims,lt,mt1,mt2]],{Res,limitdownR ,limitupR},{Ims,limitdownI,limitM},ColorFunction->"Aquamarine",AxesLabel->Automatic]];
(* Imaginary part, 1st Riemann sheet *)
Quiet[ρfig3=Plot3D[Im[ρ[Res+ I Ims,lt,mt1,mt2]],{Res,limitdownR ,limitupR},{Ims,limitM,limitupI},ColorFunction->"Aquamarine",AxesLabel->Automatic]];
Quiet[ρfig4=Plot3D[Im[ρ[Res+ I Ims,lt,mt1,mt2]],{Res,limitdownR ,limitupR},{Ims,limitdownI,limitM},ColorFunction->"Aquamarine",AxesLabel->Automatic]];

(* Here we plot *)
GraphicsRow[{Show[ρfig1,ρfig2,PlotRange->Automatic,ImageSize->Medium,PlotLabel->Style[Framed["Re ρ(s)"],FontSize->14]],Show[ρfig3,ρfig4,PlotRange->Automatic,ImageSize->Medium,PlotLabel->Style[Framed["Im ρ(s)"],FontSize->14]]},ImageSize->Large]

(* Notice the right hand cut in the right plot *)

...and ρ(s) in the complex L plane for a fixed s

scm = 2.5;

(* Ranges with the plots *)
limitdownR =0; limitupR=2*sth[mt1,mt2];limitdownI=-1;limitupI=1;limitM=0;

(* We prepare the plots *)
(* Real part, 1st Riemann sheet *)
Quiet[ρfig1=Plot3D[Re[ρ[scm,Rel+I Iml,mt1,mt2]],{Rel,limitdownR ,limitupR},{Iml,limitM,limitupI},ColorFunction->"Aquamarine",AxesLabel->Automatic,ImageSize->Medium]];
Quiet[ρfig2=Plot3D[Re[ρ[scm,Rel+I Iml,mt1,mt2]],{Rel,limitdownR ,limitupR},{Iml,limitdownI,limitM},ColorFunction->"Aquamarine",AxesLabel->Automatic,ImageSize->Medium]];
(* Imaginary part, 1st Riemann sheet *)
Quiet[ρfig3=Plot3D[Im[ρ[scm,Rel+I Iml,mt1,mt2]],{Rel,limitdownR ,limitupR},{Iml,limitM,limitupI},ColorFunction->"Aquamarine",AxesLabel->Automatic,ImageSize->Medium]];
Quiet[ρfig4=Plot3D[Im[ρ[scm,Rel+I Iml,mt1,mt2]],{Rel,limitdownR ,limitupR},{Iml,limitdownI,limitM},ColorFunction->"Aquamarine",AxesLabel->Automatic,ImageSize->Medium]];

(* Here we plot *)
GraphicsRow[{Show[ρfig1,ρfig2,PlotRange->Automatic,ImageSize->Medium,PlotLabel->Style[Framed["Re ρ(s)"],FontSize->14]],Show[ρfig3,ρfig4,PlotRange->Automatic,ImageSize->Medium,PlotLabel->Style[Framed["Im ρ(s)"],FontSize->14]]},ImageSize->Large]

(* Notice what happens when we approach Re[L]-> -1/2 *)

Single channel amplitude

Single channel reduced amplitude

(* Denominator of the amplitude in the 1st Riemann sheet*)
DENI[s_,l_,β_,α0_,α1_,m1_,m2_]:=l-α0-α1*s- β*ρ[s,l,m1,m2];

(* Reduced amplitude in the 1st Riemann sheet*)
φI[s_,l_,β_,α0_,α1_,m1_,m2_]:= β/DENI[s,l,β,α0,α1,m1,m2];

Parameters

(* Partial wave*)
l=1;

(* Masses *)
m1=0.140; m2=0.140; (* pi pi -> pi pi*)

(* Regge trajectory *) 
α0=0.5; α1=0.9; (* These parameters correspond to the ρ Regge trajectory*)

(* Residue *) 
 β=1;

(* Threshold*)
DisplayForm[{"Threshold",sth[m1,m2]}]

Reduced amplitude in the 1st Riemann sheet

(* s variable ranges, you can play with them *)
slimitdown =0;slimitup =1;

Show[GraphicsRow[
{Plot[{Re[φI[Res+ I ϵ,l,β,α0,α1,m1,m2]],Re[φI[Res-I ϵ,l,β,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{-30,30},PlotLabel->Style[Framed["Re φ(s+i0)(blue),Re φ(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m2],{Orange,Thick}}},{}}],Plot[{Im[φI[Res+ I ϵ,l,β,α0,α1,m1,m2]],Im[φI[Res-I ϵ,l,β,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{-50,50},PlotLabel->Style[Framed["Im φ(s+i0)(blue),Im φ(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m2],{Orange,Thick}}},{}}],
Plot[{Abs[φI[Res+ I ϵ,l,β,α0,α1,m1,m2]],Abs[φI[Res-I ϵ,l,β,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{0,50},PlotLabel->Style[Framed["Abs φ(s+i0)(blue),Abs φ(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m2],{Orange,Thick}}},{}}]}],ImageSize->Full]
(* The vertical orange line is the threshold *)

Pole search in the 1st Riemann sheet

DisplayForm[{"1st RS",NSolve[DENI[xs,l,β,α0,α1,m1,m2]==0,xs]}] (* only poles in the real axis are acceptable *)

Plot reduced amplitude

Quiet[φfig1=Plot3D[Re[φI[Res+ I Ims,l,β,α0,α1,m1,m2]],{Res,slimitdown,slimitup},{Ims,0,0.5},ColorFunction->"Aquamarine",PlotLabel->Style[Framed["Re φ(s)"],FontSize->14]]];
Quiet[φfig2=Plot3D[Re[φI[Res+ I Ims,l,β,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Aquamarine",PlotLabel->Style[Framed["Re φ(s)"],FontSize->14]]];
Quiet[φfig3=Plot3D[Im[φI[Res+ I Ims,l,β,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,0,0.5},ColorFunction->"Aquamarine",PlotLabel->Style[Framed["Im φ(s)"],FontSize->14]]];
Quiet[φfig4=Plot3D[Im[φI[Res+ I Ims,l,β,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Aquamarine",PlotLabel->Style[Framed["Im φ(s)"],FontSize->14]]];
Quiet[φfig5=Plot3D[Abs[φI[Res+ I Ims,l,β,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,0,0.5},ColorFunction->"Aquamarine",PlotLabel->Style[Framed["Abs φ(s)"],FontSize->14]]];
Quiet[φfig6=Plot3D[Abs[φI[Res+ I Ims,l,β,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Aquamarine"]];
GraphicsRow[{Show[φfig1,φfig2,PlotRange->Automatic,ImageSize->Medium],Show[φfig3,φfig4,PlotRange->Automatic,ImageSize->Medium],Show[φfig5,φfig6,PlotRange->Automatic,ImageSize->Medium]},ImageSize->Full]

Reduced amplitude in the 2nd Riemann sheet

(* Reduced amplitude in the 2nd Riemann sheet*)
DENII[s_,l_,β_,α0_,α1_,m1_,m2_]:=l-α0-α1*s-β*(ρ[s,l,m1,m2]+2  I Cl[s,l,m1,m2]);
φII[s_,l_,β_,α0_,α1_,m1_,m2_]:=φI[s,l,β,α0,α1,m1,m2]/(  1-2  I Cl[s,l,m1,m2] φI[s,l,β,α0,α1,m1,m2] );
φIIb[s_,l_,β_,α0_,α1_,m1_,m2_]:=β/DENII[s,l,β,α0,α1,m1,m2];
(* φII and φIIb are the same just written in two different ways, it allows to check that we are writing de denominator properly to obtain the poles (zeros) *)

(* Reduced amplitude in the first vs. in the second Riemann sheet at the physical axis *)
mt1=m1; mt2=m2;
Show[GraphicsRow[
{Plot[{Re[φI[Res+ I ϵ,l,β,α0,α1,mt1,mt2]],Re[φII[Res-I ϵ,l,β,α0,α1,mt1,mt2]]},{Res,0,slimitup},PlotRange->{-30,30},GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}},PlotStyle->{  {Blue,Thick},{Red,Dashed}},PlotLabel->Style[Framed["Re φI(s+i0) (blue), Re φII(s-i0) (red)"],FontSize->14]],Plot[{Im[φI[Res+ I ϵ,l,β,α0,α1,mt1,mt2]],Im[φII[Res-I ϵ,l,β,α0,α1,mt1,mt2]]},{Res,0,slimitup},PlotRange->{-10,50},GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}},PlotStyle->{  {Blue,Thick},{Red,Dashed}},PlotLabel->Style[Framed["Im φI(s+i0) (blue), Im φII(s-i0) (red)"],FontSize->14]],
Plot[{Abs[φI[Res+ I ϵ,l,β,α0,α1,mt1,mt2]],Abs[φII[Res-I ϵ,l,β,α0,α1,mt1,mt2]]},{Res,0,slimitup},PlotRange->{0,50},GridLines->{{{sth[mt1,mt2],{Orange,Thick}}},{}},PlotStyle->{  {Blue,Thick},{Red,Dashed}},PlotLabel->Style[Framed["Abs φI(s+i0) (blue), Abs φII(s-i0) (red)"],FontSize->14]]}],ImageSize->Full]
(* The curves entirely overlap, except below threshold, move the threshold and you will see *)

Pole search in the 2nd Riemann sheet

DisplayForm[{"2nd RS",NSolve[DENII[xs,l,β,α0,α1,m1,m2]==0,xs]}]

(* We compute the plots for the 2nd Riemann sheet ad Im s < 0*)
Quiet[φIIfig2=Plot3D[Re[φII[Res+ I Ims,l,β,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Rainbow",PlotLabel->Style[Framed["Re φ(s)"],FontSize->14]]];
Quiet[φIIfig4=Plot3D[Im[φII[Res+ I Ims,l,β,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Rainbow",PlotLabel->Style[Framed["Im φ(s)"],FontSize->14]]];
Quiet[φIIfig6=Plot3D[Abs[φII[Res+ I Ims,l,β,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Rainbow"]];
GraphicsRow[{Show[φfig1,φIIfig2,PlotRange->Automatic,ImageSize->Medium],Show[φfig3,φIIfig4,PlotRange->Automatic,ImageSize->Medium],Show[φfig5,φIIfig6,PlotRange->Automatic,ImageSize->Medium]},ImageSize->Full]

What is a Resonance?

A pole of the scattering amplitude in the unphysical Riemann sheets (s plane)

Let us see that it is also a moving pole in L and how Regge trajectories connect poles in different partial waves

Regge Trajectories

(* Our function is analytical in both s and l complex planes*)
(* Set the s limits to plot *)
slimitdown=0;
slimitup=5;
lrangedown=0; (* Do not go below L=-1/2*)
lrangeup =6; (* Depends on how much computing power you have*) 
(* Here we plot *)
Animate[Plot3D[Abs[φII[Res+I Ims,l,β,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.2,0},ColorFunction->"Rainbow",PlotLabel->Style[{"Abs φ(s)",l},FontSize->14],ImageSize->Large,Mesh->All,PlotRange->{0,20}],{l,{0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5}},AnimationRate->0.1,AnimationRunning->False]

The pole movement defines the Regge trajectory

Reggepoints[ β_, α0_, α1_, m1_, m2_,npts_]:=Table[FindRoot[ DENII[xs, l, β, α0, α1, m1, m2] == 0, {xs,sth[m1,m2]+0.1}],{l, lmin, lmax, (lmax-lmin)/(npts - 1)}];

lmin=0;lmax=6;npts=200;nptsPhys=lmax-lmin +1;
Quiet[ArrayReggePoints=xs/.Reggepoints[ β, α0, α1, m1, m2,npts]];
Quiet[ArrayReggePointsPhys=xs/.Reggepoints[ β, α0, α1, m1, m2,nptsPhys]];
ArrayReggePointsRe=Take[Re[ArrayReggePoints],{1,npts}];
ArrayReggePointsIm=Take[Im[ArrayReggePoints],{1,npts}];
ArrayReggePointsRePhys=Take[Re[ArrayReggePointsPhys],{1,nptsPhys}];
ArrayReggePointsImPhys=Take[Im[ArrayReggePointsPhys],{1,nptsPhys}];
ListPlot[{Thread[{ArrayReggePointsRe,ArrayReggePointsIm}],Thread[{ArrayReggePointsRePhys,ArrayReggePointsImPhys}]},ImageSize->Large,PlotStyle->{AbsoluteThickness[2.0], AbsolutePointSize[10]} ,AxesLabel->{Re s, Im s},PlotLabel->"Regge Trajectory"]

Notice that the Regge trajectory IS NOT LINEAR and it NEVER IS because of unitarity (as energy increases it can go linear)

ArrayLreal=Table[l,{l, lmin, lmax, (lmax-lmin)/(npts-1)}];
ArrayLPhys=Table[l,{l, lmin, lmax, (lmax-lmin)/(nptsPhys-1)}];
Quiet[ListPlot[{Thread[{ArrayReggePointsRe,ArrayLreal}],Thread[{ArrayReggePointsRePhys,ArrayLPhys}]},ImageSize->Large,PlotStyle->{AbsoluteThickness[2.0], AbsolutePointSize[10]} ,AxesLabel->{Re s=M^2, L},PlotLabel->"Chew-Frautschi Plot"]]

What can look fairly linear is the Chew-Frautschi Plot (M^2, L)

Please check PRD93, 074015 (2017), arxiv:1512.03136 for further insight

(* Let's move to the complex angular momentum plane *)
ArrayLc=Table[ArrayLreal[[i]]-β 2 I Cl[ArrayReggePoints[[i]],ArrayLreal[[i]],m1,m2],{i,1,npts}];
ArrayLcPhys=Table[ArrayLPhys[[i]]-β 2 I Cl[ArrayReggePointsPhys[[i]],ArrayLPhys[[i]],m1,m2],{i,1,nptsPhys}];ArrayLcRe =Take[Re[ArrayLc],{1,npts}];
ArrayLcIm=Take[Im[ArrayLc],{1,npts}];
ArrayLcRePhys =Take[Re[ArrayLcPhys],{1,nptsPhys}];
ArrayLcImPhys=Take[Im[ArrayLcPhys],{1,nptsPhys}];
Quiet[ListPlot[{Thread[{ArrayLcRe,ArrayLcIm}],Thread[{ArrayLcRePhys,ArrayLcImPhys}]},ImageSize->Large,PlotStyle->{AbsoluteThickness[2.0], AbsolutePointSize[10]} ,AxesLabel->{"Re L", "Im L"},PlotLabel->"Complex L plane",PlotRange->All]]

Coupled channels

Parameters

(* Partial wave*)
l=1;

(* Masses *)
m1=0.140; (* pion *)
m2=0.492; (* Kaon *)

(* Regge trajectory *) 
α0=0.5; α1=0.9;

(* Residues *)
β11=1;
β22=0.5;
β12=Sqrt[β11 β22]; (* Remember factorization *)

(* Thresholds*)
DisplayForm[{"Thresholds",sth[m1,m1],sth[m1,m2],
sth[m2,m2]}]

Reduced amplitudes in the 1st Riemann sheet

(* The Regge trajectory changes because of unitarity *)

(* Denominator of the amplitude in the 1st Riemann sheet*)
DENIcc[s_,l_,β11_,β22_,α0_,α1_,m1_,m2_]:=l-α0-α1*s- β11*ρ[s,l,m1,m1]- β22*ρ[s,l,m2,m2];

(* Reduced amplitude in the 1st Riemann sheet*)
φIcc[s_,l_,βnum_,β11_,β22_,α0_,α1_,m1_,m2_]:= βnum/DENIcc[s,l,β11,β22,α0,α1,m1,m2];

slimitdown =0;slimitup =1.2;

(* Channel 11 *)
Show[GraphicsRow[
{Plot[{Re[φIcc[Res+ I ϵ,l,β11,β11,β22,α0,α1,m1,m2]],Re[φIcc[Res-I ϵ,l,β11,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->Automatic,PlotLabel->Style[Framed["Re φ11(s+i0)(blue),Re φ11(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],Plot[{Im[φIcc[Res+ I ϵ,l,β11,β11,β22,α0,α1,m1,m2]],Im[φIcc[Res-I ϵ,l,β11,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange-> {-35,35},PlotLabel->Style[Framed["Im φ11(s+i0)(blue),Im φ11(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],
Plot[{Abs[φIcc[Res+ I ϵ,l,β11,β11,β22,α0,α1,m1,m2]],Abs[φIcc[Res-I ϵ,l,β11,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{0,40},PlotLabel->Style[Framed["Abs φ11(s+i0)(blue),Abs φ11(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}]}],ImageSize->Full]
(* Channel 12 *)
Show[GraphicsRow[
{Plot[{Re[φIcc[Res+ I ϵ,l,β12,β11,β22,α0,α1,m1,m2]],Re[φIcc[Res-I ϵ,l,β12,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->Automatic,PlotLabel->Style[Framed["Re φ12(s+i0)(blue),Re φ12(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],Plot[{Im[φIcc[Res+ I ϵ,l,β12,β11,β22,α0,α1,m1,m2]],Im[φIcc[Res-I ϵ,l,β12,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{-30,30},PlotLabel->Style[Framed["Im φ12(s+i0)(blue),Im φ12(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],
Plot[{Abs[φIcc[Res+ I ϵ,l,β12,β11,β22,α0,α1,m1,m2]],Abs[φIcc[Res-I ϵ,l,β12,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{0,30},PlotLabel->Style[Framed["Abs φ12(s+i0)(blue),Abs φ12(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}]}],ImageSize->Full]
(* Channel 22 *)
Show[GraphicsRow[
{Plot[{Re[φIcc[Res+ I ϵ,l,β22,β11,β22,α0,α1,m1,m2]],Re[φIcc[Res-I ϵ,l,β22,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->Automatic,PlotLabel->Style[Framed["Re φ22(s+i0)(blue),Re φ22(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],Plot[{Im[φIcc[Res+ I ϵ,l,β22,β11,β22,α0,α1,m1,m2]],Im[φIcc[Res-I ϵ,l,β22,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{-20,20},PlotLabel->Style[Framed["Im φ22(s+i0)(blue),Im φ22(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],
Plot[{Abs[φIcc[Res+ I ϵ,l,β22,β11,β22,α0,α1,m1,m2]],Abs[φIcc[Res-I ϵ,l,β22,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{0,20},PlotLabel->Style[Framed["Abs φ22(s+i0)(blue),Abs φ22(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}]}],ImageSize->Full]

(* The vertical orange lines are the thresholds *)

Reduced amplitude in the s complex plane

Quiet[φ11fig1=Plot3D[Re[φIcc[Res+ I Ims,l,β11,β11,β22,α0,α1,m1,m2]],{Res,slimitdown,slimitup},{Ims,0,0.5},ColorFunction->"Aquamarine",PlotLabel->Style[Framed["Re φ11(s)"],FontSize->14]]];
Quiet[φ11fig2=Plot3D[Re[φIcc[Res+ I Ims,l,β11,β11,β22,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Aquamarine",PlotLabel->Style[Framed["Re φ11(s)"],FontSize->14]]];
Quiet[φ11fig3=Plot3D[Im[φIcc[Res+ I Ims,l,β11,β11,β22,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,0,0.5},ColorFunction->"Aquamarine",PlotLabel->Style[Framed["Im φ11(s)"],FontSize->14]]];
Quiet[φ11fig4=Plot3D[Im[φIcc[Res+ I Ims,l,β11,β11,β22,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Aquamarine",PlotLabel->Style[Framed["Im φ11(s)"],FontSize->14]]];
Quiet[φ11fig5=Plot3D[Abs[φIcc[Res+ I Ims,l,β11,β11,β22,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,0,0.5},ColorFunction->"Aquamarine",PlotLabel->Style[Framed["Abs φ11(s)"],FontSize->14]]];
Quiet[φ11fig6=Plot3D[Abs[φIcc[Res+ I Ims,l,β11,β11,β22,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Aquamarine"]];
GraphicsRow[{Show[φ11fig1,φ11fig2,PlotRange->Automatic,ImageSize->Medium],Show[φ11fig3,φ11fig4,PlotRange->Automatic,ImageSize->Medium],Show[φ11fig5,φ11fig6,PlotRange->Automatic,ImageSize->Medium]},ImageSize->Full]

Poles in the 1st Riemann sheet

DisplayForm[{"1st RS",NSolve[DENIcc[xs,l,β11,β22,α0,α1,m1,m2]==0,xs]}]

Reduced amplitudes in the 2nd Riemann sheet

(* The Regge trajectory changes because of unitarity *)

(* Denominator of the amplitude in the 1st Riemann sheet*)
DENIIcc[s_,l_,β11_,β22_,α0_,α1_,m1_,m2_]:=l-α0-α1*s- β11*ρ[s,l,m1,m1]- β22*ρ[s,l,m2,m2]-2  I β11* Cl[s,l,m1,m1];

(* Reduced amplitude in the 1st Riemann sheet*)
φIIcc[s_,l_,βnum_,β11_,β22_,α0_,α1_,m1_,m2_]:= βnum/DENIIcc[s,l,β11,β22,α0,α1,m1,m2];

slimitdown =0;slimitup =1.2;

(* Channel 11 *)
Show[GraphicsRow[
{Plot[{Re[φIcc[Res+ I ϵ,l,β11,β11,β22,α0,α1,m1,m2]],Re[φIIcc[Res-I ϵ,l,β11,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->Automatic,PlotLabel->Style[Framed["Re φ11(s+i0)(blue),Re φ11^II(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],Plot[{Im[φIcc[Res+ I ϵ,l,β11,β11,β22,α0,α1,m1,m2]],Im[φIIcc[Res-I ϵ,l,β11,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange-> {-35,35},PlotLabel->Style[Framed["Im φ11(s+i0)(blue),Im φ11^II(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],
Plot[{Abs[φIcc[Res+ I ϵ,l,β11,β11,β22,α0,α1,m1,m2]],Abs[φIIcc[Res-I ϵ,l,β11,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{0,40},PlotLabel->Style[Framed["Abs φ11(s+i0)(blue),Abs φ11^II(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}]}],ImageSize->Full]
(* Channel 12 *)
Show[GraphicsRow[
{Plot[{Re[φIcc[Res+ I ϵ,l,β12,β11,β22,α0,α1,m1,m2]],Re[φIIcc[Res-I ϵ,l,β12,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->Automatic,PlotLabel->Style[Framed["Re φ12(s+i0)(blue),Re φ12^II(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],Plot[{Im[φIcc[Res+ I ϵ,l,β12,β11,β22,α0,α1,m1,m2]],Im[φIIcc[Res-I ϵ,l,β12,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{-30,30},PlotLabel->Style[Framed["Im φ12(s+i0)(blue),Im φ12^II(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],
Plot[{Abs[φIcc[Res+ I ϵ,l,β12,β11,β22,α0,α1,m1,m2]],Abs[φIIcc[Res-I ϵ,l,β12,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{0,30},PlotLabel->Style[Framed["Abs φ12(s+i0)(blue),Abs φ12^II(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}]}],ImageSize->Full]
(* Channel 22 *)
Show[GraphicsRow[
{Plot[{Re[φIcc[Res+ I ϵ,l,β22,β11,β22,α0,α1,m1,m2]],Re[φIIcc[Res-I ϵ,l,β22,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->Automatic,PlotLabel->Style[Framed["Re φ22(s+i0)(blue),Re φ22^II(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],Plot[{Im[φIcc[Res+ I ϵ,l,β22,β11,β22,α0,α1,m1,m2]],Im[φIIcc[Res-I ϵ,l,β22,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{-20,20},PlotLabel->Style[Framed["Im φ22(s+i0)(blue),Im φ22^II(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}],
Plot[{Abs[φIcc[Res+ I ϵ,l,β22,β11,β22,α0,α1,m1,m2]],Abs[φIIcc[Res-I ϵ,l,β22,β11,β22,α0,α1,m1,m2]]},{Res,slimitdown,slimitup},PlotRange->{0,20},PlotLabel->Style[Framed["Abs φ22(s+i0)(blue),Abs φ22^II(s-i0) (red)"],FontSize->14],PlotStyle->{  {Blue,Thick},{Red,Dashed}},GridLines->{{{sth[m1,m1],{Orange,Thick}},{sth[m2,m2],{Orange,Thick}}},{}}]}],ImageSize->Full]

(* The vertical orange lines are the thresholds *)
(* You see how 1st and 2nd RS beautifully match at the real axis, except below and above threshold, move the threshold and you will see *)

Amplitudes and pole search in the 2nd Riemann sheet

DisplayForm[{"2nd RS",NSolve[DENIIcc[xs,l,β11,β22,α0,α1,m1,m2]==0,xs]}]

(* We compute the plots for the 2nd Riemann sheet ad Im s < 0*)
Quiet[φ11IIfig2=Plot3D[Re[φIIcc[Res+ I Ims,l,β11,β11,β22,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Rainbow",PlotLabel->Style[Framed["Re φ(s)"],FontSize->14]]];
Quiet[φ11IIfig4=Plot3D[Im[φIIcc[Res+ I Ims,l,β11,β11,β22,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Rainbow",PlotLabel->Style[Framed["Im φ(s)"],FontSize->14]]];
Quiet[φ11IIfig6=Plot3D[Abs[φIIcc[Res+ I Ims,l,β11,β11,β22,α0,α1,m1,m2]],{Res,slimitdown,slimitup },{Ims,-0.5,0},ColorFunction->"Rainbow"]];
GraphicsRow[{Show[φ11fig1,φ11IIfig2,PlotRange->Automatic,ImageSize->Medium],Show[φ11fig3,φ11IIfig4,PlotRange->Automatic,ImageSize->Medium],Show[φ11fig5,φ11IIfig6,PlotRange->Automatic,ImageSize->Medium]},ImageSize->Full]

Regge Trajectories for the coupled channel case

ReggepointsCC[ β11_,β22_, α0_, α1_, m1_, m2_,npts_]:=Table[FindRoot[ DENIIcc[xs, l, β11,β22, α0, α1, m1, m2] == 0, {xs,sth[m1,m2]+0.1}],{l, lmin, lmax, (lmax-lmin)/(npts - 1)}];

lminCC=0;lmaxCC=6;nptsCC=200;nptsPhysCC=lmax-lmin +1;
Quiet[ArrayReggePointsCC=xs/.ReggepointsCC[ β11,β22, α0, α1, m1, m2,nptsCC]];
Quiet[ArrayReggePointsPhysCC=xs/.ReggepointsCC[ β11,β22, α0, α1, m1, m2,nptsPhysCC]];
ArrayReggePointsReCC=Take[Re[ArrayReggePointsCC],{1,nptsCC}];
ArrayReggePointsImCC=Take[Im[ArrayReggePointsCC],{1,nptsCC}];
ArrayReggePointsRePhysCC=Take[Re[ArrayReggePointsPhysCC],{1,nptsPhysCC}];
ArrayReggePointsImPhysCC=Take[Im[ArrayReggePointsPhysCC],{1,nptsPhysCC}];
ListPlot[{Labeled[Thread[{ArrayReggePointsRe,ArrayReggePointsIm}],"Single channel"],Thread[{ArrayReggePointsRePhys,ArrayReggePointsImPhys}],Labeled[Thread[{ArrayReggePointsReCC,ArrayReggePointsImCC}],"Coupled channel"],Thread[{ArrayReggePointsRePhysCC,ArrayReggePointsImPhysCC}]},ImageSize->Large,PlotStyle->{AbsoluteThickness[2.0], AbsolutePointSize[10]} ,AxesLabel->{Re s, Im s},PlotLabel->"Regge Trajectory",PlotRange->All]

ArrayLrealCC=Table[l,{l, lminCC, lmaxCC, (lmaxCC-lminCC)/(nptsCC-1)}];
ArrayLPhysCC=Table[l,{l, lminCC, lmaxCC, (lmaxCC-lminCC)/(nptsPhysCC-1)}];
Quiet[ListPlot[{Labeled[Thread[{ArrayReggePointsRe,ArrayLreal}],"Single channel"],Thread[{ArrayReggePointsRePhys,ArrayLPhys}],Labeled[Thread[{ArrayReggePointsReCC,ArrayLrealCC}],"Coupled channel"],Thread[{ArrayReggePointsRePhysCC,ArrayLPhysCC}]},ImageSize->Large,PlotStyle->{AbsoluteThickness[2.0], AbsolutePointSize[10]} ,AxesLabel->{Re s=M^2, L},PlotLabel->"Chew-Frautschi Plot"]]

Notice how similar are the Chew-Frautschi plots... and how different are the Regge trajectories!!

(* Let's move to the complex angular momentum plane *)

ArrayLcCC=Table[ArrayLrealCC[[i]]-β 2 I Cl[ArrayReggePointsCC[[i]],ArrayLrealCC[[i]],m1,m2],{i,1,nptsCC}];
ArrayLcPhysCC=Table[ArrayLPhysCC[[i]]-β 2 I Cl[ArrayReggePointsPhysCC[[i]],ArrayLPhysCC[[i]],m1,m2],{i,1,nptsPhysCC}];ArrayLcReCC =Take[Re[ArrayLcCC],{1,nptsCC}];
ArrayLcImCC=Take[Im[ArrayLcCC],{1,nptsCC}];
ArrayLcRePhysCC =Take[Re[ArrayLcPhysCC],{1,nptsPhysCC}];
ArrayLcImPhysCC=Take[Im[ArrayLcPhysCC],{1,nptsPhysCC}];
Quiet[ListPlot[{Labeled[Thread[{ArrayLcRe,ArrayLcIm}],"Single channel"],Thread[{ArrayLcRePhys,ArrayLcImPhys}],Labeled[Thread[{ArrayLcReCC,ArrayLcImCC}],"Coupled channel"],Thread[{ArrayLcRePhysCC,ArrayLcImPhysCC}]},ImageSize->Large,PlotStyle->{AbsoluteThickness[2.0], AbsolutePointSize[10]} ,AxesLabel->{"Re L", "Im L"},PlotLabel->"Complex L plane"]]