Construction of the 2 -> 2 Scattering Matrices:

In the following we describe the scalar potential for the LRSM. Later we construct the scattering matrices q-charged 2-particle states by properly specifying the basis.
(Uncomment and run the command “Eigenvalues[X]” to obtain the eigenvalues of the matrix ‘X’.)

Defining the scalar potential

Clear["Global`*"]

ϕ=(({{(1/(Sqrt[2]))(ϕ10r+I  ϕ10i), ϕ1p}, {ϕ2m, (1/(Sqrt[2]))(ϕ20r+I ϕ20i)}}));ΔR=(({{((δRp)/(Sqrt[2])), δRpp}, {(1/(Sqrt[2]))(δR0r+I δR0i), -((δRp)/(Sqrt[2]))}}));ΔL=(({{((δLp)/(Sqrt[2])), δLpp}, {(1/(Sqrt[2]))(δL0r+I δL0i), -((δLp)/(Sqrt[2]))}}));
σ=(({{0, -I}, {I, 0}}));ϕt = σ.ϕ.σ;
V4 = (Subscript[λ, 1]) (((Tr[ϕ\[ConjugateTranspose].ϕ]))^2)+(Subscript[λ, 2])((((Tr[ϕt.ϕ\[ConjugateTranspose]]))^2)+(((Tr[ϕt\[ConjugateTranspose].ϕ]))^2))+(Subscript[λ, 3])((Tr[ϕt.ϕ\[ConjugateTranspose]])(Tr[ϕt\[ConjugateTranspose].ϕ]))+(Subscript[λ, 4])((Tr[ϕ\[ConjugateTranspose].ϕ])((Tr[ϕt\[ConjugateTranspose].ϕ])+(Tr[ϕt.ϕ\[ConjugateTranspose]])))+(Subscript[λ, 5])((((Tr[ΔL.ΔL\[ConjugateTranspose]]))^2)+(((Tr[ΔR.ΔR\[ConjugateTranspose]]))^2))+(Subscript[λ, 6])( (Tr[ΔL.ΔL])(Tr[ΔL\[ConjugateTranspose].ΔL\[ConjugateTranspose]])+ (Tr[ΔR.ΔR])(Tr[ΔR\[ConjugateTranspose].ΔR\[ConjugateTranspose]]))+(Subscript[λ, 7]) (Tr[ΔL.ΔL\[ConjugateTranspose]]Tr[ΔR.ΔR\[ConjugateTranspose]])+(Subscript[λ, 8])(Tr[ΔL.ΔL]Tr[ΔR\[ConjugateTranspose].ΔR\[ConjugateTranspose]]+ Tr[ΔL\[ConjugateTranspose].ΔL\[ConjugateTranspose]]Tr[ΔR.ΔR])+(Subscript[λ, 9]) Tr[ϕ\[ConjugateTranspose].ϕ](Tr[ΔL\[ConjugateTranspose].ΔL\[ConjugateTranspose]]+Tr[ΔR\[ConjugateTranspose].ΔR\[ConjugateTranspose]])+((Subscript[λ, 10])+I (Subscript[λ, 11]))((Tr[ϕt\[ConjugateTranspose].ϕ])(Tr[ΔR.ΔR\[ConjugateTranspose]])+(Tr[ϕ\[ConjugateTranspose].ϕt])(Tr[ΔL.ΔL\[ConjugateTranspose]]))+((Subscript[λ, 10])-I (Subscript[λ, 11]))((Tr[ϕt\[ConjugateTranspose].ϕ])(Tr[ΔL.ΔL\[ConjugateTranspose]])+(Tr[ϕt.ϕ\[ConjugateTranspose]])(Tr[ΔR.ΔR\[ConjugateTranspose]]))+(Subscript[λ, 12])(Tr[ϕ\[ConjugateTranspose].ϕ.ΔL.ΔL\[ConjugateTranspose]]+Tr[ϕ\[ConjugateTranspose].ϕ.ΔR.ΔR\[ConjugateTranspose]]);

replace={Conjugate[ϕ10r]-> ϕ10r,Conjugate[ϕ10i]-> ϕ10i,Conjugate[ϕ20r]-> ϕ20r,Conjugate[ϕ20i]-> ϕ20i,Conjugate[ϕ1p]-> ϕ1m,Conjugate[ϕ2m]-> ϕ2p,Conjugate[δR0r]-> δR0r,Conjugate[δR0i]-> δR0i,Conjugate[δRp]-> δRm,Conjugate[δRpp]-> δRmm,Conjugate[δL0r]-> δL0r,Conjugate[δL0i]-> δL0i,Conjugate[δLp]-> δLm,Conjugate[δLpp]-> δLmm,Conjugate[I ϕ10i+ϕ10r]-> {- I ϕ10i+ϕ10r},Conjugate[I ϕ20i+ϕ20r]-> {- I ϕ20i+ϕ20r}};

(*Potential in terms of the fields *)
Vtot=Simplify[V4/.replace];

Construction of scattering amplitude matrix

For neutral two-particle states

NCColumn={
ϕ1p ϕ1m,ϕ1p ϕ2m,ϕ1p δRm ,ϕ1p δLm (*4,4*),
ϕ2p ϕ1m, ϕ2p ϕ2m,ϕ2p δRm , ϕ2p  δLm (*4,8*),
δRp ϕ1m,δRp ϕ2m,δRp δRm,δRp δLm(*4,12*),
δLp ϕ1m,δLp ϕ2m,δLp δRm,δLp δLm(*4,16*),
ϕ10r ϕ10r,ϕ10r ϕ10i ,ϕ10r ϕ20r ,ϕ10r ϕ20i ,ϕ10r  δR0r,ϕ10r  δR0i ,ϕ10r  δL0r,ϕ10r  δL0i  (*8,24*),
ϕ10i ϕ10i,ϕ10i ϕ20r ,ϕ10i ϕ20i ,ϕ10i δR0r ,ϕ10i δR0i ,ϕ10i δL0r ,ϕ10i  δL0i(*7,31*), 
ϕ20r  ϕ20r,ϕ20r  ϕ20i,ϕ20r  δR0r,ϕ20r δR0i ,ϕ20r  δL0r,ϕ20r δL0i  (*6,37*),
ϕ20i  ϕ20i,ϕ20i δR0r ,ϕ20i δR0i ,ϕ20i δL0r  ,ϕ20i  δL0i(*5,42*), 
δR0r  δR0r,δR0r δR0i ,δR0r  δL0r,δR0r δL0i (*4,46*),
δR0i δR0i ,δR0i δL0r  ,δR0i  δL0i (*3,49*),
δL0r  δL0r,δL0r  δL0i (*2,51*),
δL0i δL0i(*1,52*), 
δRpp δRmm ,δRpp δLmm (*2,54*),
δLpp  δRmm,δLpp δLmm (*2,56*)
};
NCRow={
ϕ1m ϕ1p,ϕ1m ϕ2p,ϕ1m δRp ,ϕ1m δLp (*4,4*),
ϕ2m ϕ1p, ϕ2m ϕ2p,ϕ2m δRp , ϕ2m  δLp (*4,8*),
δRm ϕ1p,δRm ϕ2p,δRm δRp,δRm δLp(*4,12*),
δLm ϕ1p,δLm ϕ2p,δLm δRp,δLm δLp(*4,16*),
ϕ10r ϕ10r,ϕ10r ϕ10i ,ϕ10r ϕ20r ,ϕ10r ϕ20i ,ϕ10r  δR0r,ϕ10r  δR0i ,ϕ10r  δL0r,ϕ10r  δL0i  (*8,24*),
ϕ10i ϕ10i,ϕ10i ϕ20r ,ϕ10i ϕ20i ,ϕ10i δR0r ,ϕ10i δR0i ,ϕ10i δL0r ,ϕ10i  δL0i(*7,31*),
ϕ20r  ϕ20r,ϕ20r  ϕ20i,ϕ20r  δR0r,ϕ20r δR0i ,ϕ20r  δL0r,ϕ20r δL0i  (*6,37*),
ϕ20i  ϕ20i,ϕ20i δR0r ,ϕ20i δR0i ,ϕ20i δL0r  ,ϕ20i  δL0i(*5,42*),
δR0r  δR0r,δR0r δR0i ,δR0r  δL0r,δR0r δL0i (*4,46*),
δR0i δR0i ,δR0i δL0r  ,δR0i  δL0i (*3,49*),
δL0r  δL0r,δL0r  δL0i (*2,51*),
δL0i δL0i(*1,52*), 
δRmm δRpp ,δRmm δLpp (*2,54*),
δLmm  δRpp,δLmm δLpp (*2,56*)
};

The symmetry factor is introduced to avoid double counting:

NCfacList[i_,j_]:=FactorList[NCColumn[[i]] NCRow[[j]]];
NCSymFac[i_,j_]:=(If[Length[NCfacList[i,j]]==2,Return[Factorial[4]]];
If[Length[NCfacList[i,j]]==3,If[(NCfacList[i,j])[[2,2]]==2&&(NCfacList[i,j])[[3,2]]==2,Return[4],Return[Factorial[3]]]];
If[Length[NCfacList[i,j]]==4,Return[2]];
If[Length[NCfacList[i,j]]==5,Return[1]]);

NCcoeff[i_,j_]:=Coefficient[Vtot,NCColumn[[i]] NCRow[[j]]];
NCMatFunc[i_,j_]:=NCSymFac[i,j] NCcoeff[i,j];

For[i=1;NCklist={},i≤ 56,i++,For[j=1;NCjlist={},j≤ 56,j++,AppendTo[NCjlist,{NCMatFunc[i,j]}]];AppendTo[NCklist,Flatten[NCjlist]]]
(*NCklist//MatrixForm*)
(*Eigenvalues[NCklist]*)

For singly charged two-particle states

CCColumn={
ϕ1p ϕ10r ,ϕ1p ϕ10i,ϕ1p ϕ20r,ϕ1p ϕ20i ,ϕ1p δR0r,ϕ1p δR0i,ϕ1p δL0r,ϕ1p δL0i(*8,8*),
ϕ2p  ϕ10r,ϕ2p  ϕ10i,ϕ2p  ϕ20r,ϕ2p ϕ20i ,ϕ2p δR0r ,ϕ2p  δR0i,ϕ2p  δL0r,ϕ2p  δL0i (*8,16*), 
δRp ϕ10r , δRp  ϕ10i, δRp  ϕ20r, δRp  ϕ20i, δRp  δL0r,δRp δL0i  ,δRp δR0r , δRp  δR0i(*8,24*),
δLp ϕ10r,δLp ϕ10i,δLp ϕ20r ,δLp ϕ20i,δLp δR0r  ,δLp δR0i,δLp δL0r,δLp δL0i (*8,32*),
ϕ1p δRmm,ϕ1p δLmm(*2,34*), 
ϕ2p δRmm,ϕ2p δLmm(*2,36*),
δRp δRmm,δRp δLmm(*2,38*),
δLp δRmm,δLp δLmm(*2,40*)
};
CCRow={ϕ1m ϕ10r ,ϕ1m ϕ10i,ϕ1m ϕ20r,ϕ1m ϕ20i ,ϕ1m δR0r,ϕ1m δR0i,ϕ1m δL0r,ϕ1m δL0i(**),
ϕ2m  ϕ10r,ϕ2m  ϕ10i,ϕ2m  ϕ20r,ϕ2m ϕ20i ,ϕ2m δR0r ,ϕ2m  δR0i,ϕ2m  δL0r,ϕ2m  δL0i(**), 
δRm ϕ10r , δRm  ϕ10i, δRm  ϕ20r, δRm  ϕ20i, δRm  δL0r,δRm δL0i, δRm δR0r , δRm  δR0i(**),
δLm ϕ10r,δLm ϕ10i,δLm ϕ20r ,δLm ϕ20i,δLm δR0r  ,δLm δR0i ,δLm δL0r,δLm δL0i (**),
ϕ1m δRpp,ϕ1m δLpp (**),
 ϕ2m δRpp,ϕ2m δLpp (**),
δRm δRpp,δRm δLpp (**),
δLm δRpp,δLm δLpp (**)
};

Symmetry factor:

CCfacList[i_,j_]:=FactorList[CCColumn[[i]] CCRow[[j]]];
CCSymFac[i_,j_]:=(If[Length[CCfacList[i,j]]==2,Return[Factorial[4]]];
If[Length[CCfacList[i,j]]==3,If[(CCfacList[i,j])[[2,2]]==2&&(CCfacList[i,j])[[3,2]]==2,Return[4],Return[Factorial[3]]]];
If[Length[CCfacList[i,j]]==4,Return[2]];
If[Length[CCfacList[i,j]]==5,Return[1]]);

CCcoeff[i_,j_]:=Coefficient[Vtot,CCColumn[[i]] CCRow[[j]]];
CCMatFunc[i_,j_]:=CCSymFac[i,j] CCcoeff[i,j];

For[i=1;CCklist={},i≤ 40,i++,For[j=1;CCjlist={},j≤ 40,j++,AppendTo[CCjlist,{CCMatFunc[i,j]}]]
AppendTo[CCklist,Flatten[CCjlist]]]
(*CCklist//MatrixForm;*)

(*Eigenvalues[CCklist]*)

For doubly charged two-particle states

(*Eigenvalues[NCklist]*)

DCColumn={
ϕ1p ϕ1p, ϕ1p  ϕ2p, ϕ1p δRp ,ϕ1p δLp (*4,4*),
ϕ2p  ϕ2p, ϕ2p δRp ,ϕ2p δLp (*3,7*), 
δRp δRp ,δRp δLp (*2,9*),
δLp δLp (*1,10*), 
δRpp ϕ10r, δRpp ϕ10i,δRpp ϕ20r,δRpp ϕ20i,δRpp δR0r ,δRpp δR0i,δRpp δL0r,δRpp δL0i (*8,18*), 
δLpp ϕ10r, δLpp ϕ10i,δLpp ϕ20r,δLpp ϕ20i,δLpp δR0r ,δLpp δR0i,δLpp δL0r,δLpp δL0i(*8,26*) 
};
DCRow={
ϕ1m ϕ1m, ϕ1m  ϕ2m, ϕ1m δRm ,ϕ1m δLm (**),
ϕ2m  ϕ2m, ϕ2m δRm ,ϕ2m δLm (**),
δRm δRm ,δRm δLm (**), 
δLm δLm (**),
δRmm ϕ10r, δRmm ϕ10i,δRmm ϕ20r,δRmm ϕ20i,δRmm δR0r ,δRmm δR0i,δRmm δL0r,δRmm δL0i (**),
δLmm ϕ10r, δLmm ϕ10i,δLmm ϕ20r,δLmm ϕ20i,δLmm δR0r ,δLmm δR0i,δLmm δL0r,δLmm δL0i(**)
 };

Symmetry factor:

DCfacList[i_,j_]:=FactorList[DCColumn[[i]] DCRow[[j]]];
DCSymFac[i_,j_]:=(If[Length[DCfacList[i,j]]==2,Return[Factorial[4]]];
If[Length[DCfacList[i,j]]==3,If[(DCfacList[i,j])[[2,2]]==2&&(DCfacList[i,j])[[3,2]]==2,Return[4],Return[Factorial[3]]]];
If[Length[DCfacList[i,j]]==4,Return[2]];
If[Length[DCfacList[i,j]]==5,Return[1]]);

DCcoeff[i_,j_]:=Coefficient[Vtot,DCColumn[[i]] DCRow[[j]]];
DCMatFunc[i_,j_]:=DCSymFac[i,j] DCcoeff[i,j];

For[i=1;DCklist={},i≤ 26,i++,For[j=1;DCjlist={},j≤ 26,j++,AppendTo[DCjlist,{DCMatFunc[i,j]}]];AppendTo[DCklist,Flatten[DCjlist]]]
(*DCklist//MatrixForm;*)

Eigenvalues[DCklist]

For triply charged two-particle states

TCColumn={ ϕ1p  δRpp,ϕ1p δLpp ,ϕ2p δRpp ,ϕ2p δLpp,δRp  δRpp,δRp  δLpp,δLp  δRpp,δLp  δLpp};
TCRow={ ϕ1m  δRmm,ϕ1m δLmm ,ϕ2m δRmm ,ϕ2m δLmm,δRm  δRmm,δRm  δLmm,δLm  δRmm,δLm  δLmm};

Symmetry factor:

TCfacList[i_,j_]:=FactorList[TCColumn[[i]] TCRow[[j]]];
TCSymFac[i_,j_]:=(If[Length[TCfacList[i,j]]==2,Return[Factorial[4]]];
If[Length[TCfacList[i,j]]==3,If[(TCfacList[i,j])[[2,2]]==2&&(TCfacList[i,j])[[3,2]]==2,Return[4],Return[Factorial[3]]]];
If[Length[TCfacList[i,j]]==4,Return[2]];
If[Length[TCfacList[i,j]]==5,Return[1]]);

TCcoeff[i_,j_]:=Coefficient[Vtot,TCColumn[[i]] TCRow[[j]]];
TCMatFunc[i_,j_]:=TCSymFac[i,j] TCcoeff[i,j];

For[i=1;TCklist={},i≤8 ,i++,For[j=1;TCjlist={},j≤ 8,j++,AppendTo[TCjlist,{TCMatFunc[i,j]}]];AppendTo[TCklist,Flatten[TCjlist]]]
(*TCklist//MatrixForm*)
Eigenvalues[TCklist]

For quartically charged two-particle states

QCColumn={ δRpp δRpp,δRpp δLpp ,δLpp δLpp };
QCRow={δRmm δRmm,δRmm δLmm ,δLmm δLmm};

Symmetry factor:

QCfacList[i_,j_]:=FactorList[QCColumn[[i]] QCRow[[j]]];
QCSymFac[i_,j_]:=(If[Length[QCfacList[i,j]]==2,Return[Factorial[4]]];
If[Length[QCfacList[i,j]]==3,If[(QCfacList[i,j])[[2,2]]==2&&(QCfacList[i,j])[[3,2]]==2,Return[4],Return[Factorial[3]]]];
If[Length[QCfacList[i,j]]==4,Return[2]];
If[Length[QCfacList[i,j]]==5,Return[1]]);
QCcoeff[i_,j_]:=Coefficient[Vtot,QCColumn[[i]] QCRow[[j]]];
QCMatFunc[i_,j_]:=QCSymFac[i,j] QCcoeff[i,j];

For[i=1;QCklist={},i≤3 ,i++,For[j=1;QCjlist={},j≤ 3,j++,AppendTo[QCjlist,{QCMatFunc[i,j]}]];AppendTo[QCklist,Flatten[QCjlist]]]
(*QCklist//MatrixForm*)
Eigenvalues[QCklist]