Tree-Level BCFW Recursion in Momentum-Twistors

Goal: define a function “rAmp[n_,k_]” which returns a List of products of R-invariants whose sum represents the n-point (N^k)MHV amplitude

Recusion as Replacement

Strategy: define a replacement rule, which acts on a (forever!!!) undefined symbol “amp[k_][legList_]”, and repeatedly replace this with the RHS of the BCFW formula.

bcfwRecurse=amp[k_][legList_]:>Which[
(*MHV amplitudes are 1*)
k==0,1,
(*worse than MHV-bar amplitudes are 0*)
k>Length[legList]-4,0,
(*ELSE*)
True,
amp[k][legList[[1;;-2]]]+Sum[bcfwBridge[bridgeIndices][legList],{bridgeIndices,bcfwPartitions[Length[legList],k]}]
];

Partitions across the bridge

Goal: for any given {n,k}, we want a list of “termIndices” that appear in the summand—formatted according to:
{{nL,kL},{nR,kR}}.
Recall that
n=nL+nR-2;
k=kL+kR+1.

With[{n=8,k=2},
Select[Join@@Table[{{n-nR+2,k-kR-1},{nR,kR}},{nR,4,n-1},{kR,0,nR-4}],(And@@nonVanishingQ@@@#)&]
]
(And@@nonVanishingQ@@@#)&/@%

Let’s cull all the vanishing terms from the list

nonVanishingQ[n_,k_]:=Or[0≤k≤n-4,{n,k}=={3,0}];

bcfwPartitions[n_,k_]:=Select[Join@@Table[{{n-nR+2,k-kR-1},{nR,kR}},{nR,4,n-1},{kR,0,nR-4}],(And@@nonVanishingQ@@@#)&]

Application: count the terms in BCFW

With[{n=8,k=2},
termsInBCFW[n-1,k]+Total[Times@@termsInBCFW@@@#&/@bcfwPartitions[n,k]]
]

termsInBCFW[n_,k_]:=Which[k==0,1,k>n-4,0,True,termsInBCFW[n-1,k]+Total[Times@@termsInBCFW@@@#&/@bcfwPartitions[n,k]]];

termsInBCFW[8,2]

termsInBCFW[14,5]//Timing

Speed things up by saving intermediate results to memory

ClearAll[termsInBCFW];
termsInBCFW[n_,k_]:=(termsInBCFW[n,k]=Which[k==0,1,k>n-4,0,True,termsInBCFW[n-1,k]+Total[Times@@termsInBCFW@@@#&/@bcfwPartitions[n,k]]]);

termsInBCFW[14,5]//Timing

termsInBCFW[14,5]//Timing

Building the BCFW Bridge

We had a function called “bcfwBridge[termIndices_][legList_]”; we’d like this to return the product of lower-point amplitudes (and the R-invariant!) that appears in the summand of the recursion formula.

bcfwBridge[{{nL_,kL_},{nR_,kR_}}][legList_]:=
Block[{n=Length[legList],nHat,a=nL,aHat,leftLegs=legList[[1;;nL]],rightLegs=legList[[nL-1;;-1]]},
aHat=cap[legList[[{a,a-1}]],legList[[{n-1,n,1}]]];
nHat=cap[legList[[{n,n-1}]],legList[[{a-1,a,1}]]];
leftLegs=ReplacePart[leftLegs,{-1->aHat}];
rightLegs=ReplacePart[rightLegs,{1->aHat,-1->nHat}];
Times@@{amp[kL][leftLegs],R@@legList[[{1,a-1,a,n-1,n}]],amp[kR][rightLegs]}
]

Recall that 
aHat:=(a,a-1)⋂(n-1,n,1):=cap[{a,a-1},{n-1,n,1}]

With[{n=8,k=2},
Table[bcfwBridge[termIndices][Range[n]],{termIndices,bcfwPartitions[n,k]}]
]

Putting everything together

We are now ready to define the function “rAmp[n_,k_]” (!)

With[{n=7,k=2},
MonomialList[amp[k][Range[n]]//.bcfwRecurse]
]

Remember: I wanted a List of terms (aka “Monomials”)

rAmp[n_,k_]:=MonomialList[amp[k][Range[n]]//.bcfwRecurse];

Mixing it up

Let’s randomly choose which pair of legs to use for {1,n} by randomly rotating every amp[][] before it gets recursed

With[{legList=Range[7]},
RotateLeft[legList,RandomInteger[{0,Length[legList]}]]
]

randomRotateAmp=amp[k_][legList_]:>amp[k][RotateLeft[legList,RandomInteger[{0,Length[legList]}]]];

bcfwRecurseRandom=amp[k_][legList_]:>((amp[k][legList]/.randomRotateAmp)/.bcfwRecurse);

rAmpRandom[n_,k_]:=MonomialList[amp[k][Range[n]]//.bcfwRecurseRandom];

rAmpRandom[8,2];

Teaching Mathematica Supersymmetry (cleverly)

Goal: take the products or R-invariants and return “standardized” superfunctions:
f×δ(C·η)↔{f,C}↔{bosonicBit,etaCoeffs}

Bosonic bit (trivial)

R[a,b,c,d,e]/.R[x__]:>1/Times@@ab@@@Partition[List[x],4,1,1]

bosonicRule=R[x__]:>1/Times@@ab@@@Partition[List[x],4,1,1];

Eta coefficients (less trivial)

rAmp[8,2]//Last//Last
%//.{R[x__]:>((ab@@@RotateLeft[Partition[{x},4,1,1]]).η/@{x}),η[cap[{a_,b_},{c__}]]:>η[a]ab[b,c]+η[b]ab[c,a]}
Normal[CoefficientArrays[%,η/@Range[8]][[2]]];

Recall that 
(a,b)⋂(c,d,e):=z[a]⟨b,c,d,e⟩+z[b]⟨c,d,e,a⟩

We want the matrix of “Coefficients” of the η’s

?CoefficientArrays

Putting everything together

(*Warning: "exprn" must be a product of Rs*)
toSuperfunction[n_][exprn_]:=
Block[{bosonicBit=(exprn//.bosonicRule),etaCoeffs,rList=Cases[exprn,R[__],{0,∞}]},
etaCoeffs=Normal[CoefficientArrays[#,η/@Range[n]][[2]]]&/@(rList//.{R[x__]:>((ab@@@RotateLeft[Partition[{x},4,1,1]]).η/@{x}),η[cap[{a_,b_},{c__}]]:>η[a]ab[b,c]+η[b]ab[c,a]});
{bosonicBit,etaCoeffs}
]

superAmp[n_,k_]:=toSuperfunction[n]/@rAmp[n,k];
superAmpRandom[n_,k_]:=toSuperfunction[n]/@rAmpRandom[n,k];

superAmp[6,2]//First
Last[%]//MatrixForm

Extracting components

We want to define a reference component function to test our code. 
We’re going to take the “split” (helicity) component function—which corresponds to taking the first k×k minor four times

splitComponent[{bosonicBit_,etaCoeffs_}]:=
Block[{k=Length[etaCoeffs],minorList},
minorList=Range[k]&/@Range[4];
bosonicBit*Product[Det[etaCoeffs[[All,minor]]],{minor,minorList}]
]

splitComponent/@superAmpRandom[8,1]

Evaluation and Verification(!)

The kinematic data is specified by a (4×n) matrix of momentum twistors.
Note: an unconstrained set of momentum twistors always defines momentum-conserving, on-shell (massless) four-momenta. This means that may as well define more than we may need.

momentumTwistorData=RandomInteger[{1,50},{20,4}]

Teach (very cautiously) Mathematica what “ab[__]” means

evaluate[exprn_]:=Block[{out},
ab[x___Integer]:=ab[x]=Det[momentumTwistorData[[{x}]]];
ab[x___,cap[{a_,b_},{c__}],y___]:=ab[x,a,y]ab[b,c]+ab[x,b,y]ab[c,a];
out=exprn;
ClearAll[ab];
out
]

evaluate[splitComponent/@superAmp[6,2]]

splitComponent/@superAmpRandom[6,2]
evaluate[%]

With[{n=10,k=5},
evaluate[splitComponent/@superAmpRandom[n,k]];
]
Total[%]

Cancellation of spurious poles

Recall that the only physical poles of a planar amplitude are Mandelstams of the form s[a,..b-1]∝⟨a-1,a,b-1,b⟩

With[{n=8},
Times@@ab@@@DeleteCases[Join@@@Subsets[Partition[Range[n],2,1,1],{2}],
{x___,y_,z___,y_,w___}]
]

physicalPoles[n_]:=Times@@ab@@@DeleteCases[Join@@@Subsets[Partition[Range[n],2,1,1],{2}],
{x___,y_,z___,y_,w___}]

With[{n=10,k=3},
evaluate[physicalPoles[n]]*evaluate[splitComponent/@superAmpRandom[n,k]];
]
IntegerQ/@%
IntegerQ[Total[%%]]

Homework

Improve efficiency of everything above

Generalize the component extraction

Uplift this to momentum space