Tree-Level Recursion in Momentum-Twistor Space

Goal: define a function rAmp[n_,k_] which gives the n-point (N^k)MHV amplitude, 𝒜(1,...,n), given as a List of products of R-invariants. 

Important: verify that we are right.

Recursion via replacement

Strategy: because the recursion relations involve amplitudes with miscellaneous arguments including geometrically defined twistors, it’s a good idea to carry this information explicitly. 

So consider a symbol amp[k_][legLabels_], and then use Replace to repeatedly implement recursion

bcfwRecurse=amp[k_][legLabels_]:>Which[
(*Parke-Taylor amplitudes are 1*)
k==0,1,
(*worse-than PT-bar are 0*)
k>Length[legLabels]-4,0,
(*ELSE*)
True,
amp[k][legLabels[[1;;-2]]]+
Sum[bcfwBridge[termIndex][legLabels],{termIndex,bcfwPartitions[Length[legLabels],k]}]
];

Nota bene: we need to define bcfwBridge[][] and bcfwPartitions.

Partitions across the bridge

Goal: given {n,k}, we want a list of terms labeling the pieces in the summand.
Let’s format this data according to: {{nL,kL},{nR,kR}}.

Recall:
nL+nR-2=n;
kL+kR+1=k;

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@@@#&]
]
Last[%]
And@@nonVanishingQ@@@%

Let’s cull from this list all the terms that involve vanishing amplitudes.

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: counting terms in the recursion

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[14,5]//Timing

Note: we can speed things up a lot by storing intermediate values to memory using SetDelayed[Set[]]

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

Given the term indices output by bcfwPartitions, we’d like to give the term in the summand of the recursion relations.

bcfwBridge[{{nL_,kL_},{nR_,kR_}}][legLabels_]:=
Block[{n=Length[legLabels],nHat,a=nL,aHat,leftLabels,rightLabels},
aHat=cap[legLabels[[{a-1,a}]],legLabels[[{n-1,n,1}]]];
nHat=cap[legLabels[[{n-1,n}]],legLabels[[{a-1,a,1}]]];
leftLabels=ReplacePart[legLabels[[1;;nL]],{-1->aHat}];
rightLabels=ReplacePart[legLabels[[nL-1;;-1]],{1->aHat,-1->nHat}];
Times@@{amp[kL][leftLabels],R@@legLabels[[{1,a-1,a,n-1,n}]],amp[kR][rightLabels]}
];

Recall that 
aHat:=(a-1 a)⋂(n-1 n 1):=”cap[{a-1,a},{n-1,n,1}]” and similarly for nHat.
(Recall that the LATEX command for “⋂” is “\cap”.)

With[{n=7,k=2},
Table[bcfwBridge[termIndex][Range[7]],{termIndex,bcfwPartitions[n,k]}]
]

Putting things together

Recall that we wanted rAmp to return a List of terms (aka Monomials)

With[{n=8,k=1},
MonomialList[amp[k][Range[n]]//.bcfwRecurse]
]

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

Mixing things up a bit

If we did things correctly, all amplitudes should be cyclically invariant. Let’s randomize this.

With[{legLabels=Range[7]},
RotateLeft[legLabels,RandomInteger[Length[legLabels]-1]]
]

rotateAmpRandom=amp[k_][legLabels_]:>amp[k][RotateLeft[legLabels,RandomInteger[Length[legLabels]-1]]];

bcfwRecurseRandom=amp[k_][legLabels_]:>((amp[k][legLabels]/.rotateAmpRandom)/.bcfwRecurse);

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

rAmpRandom[8,2]

Teaching Mathematica Supersymmetry (non-stupidly)

Standard superfunctions

Goal: convert the products of R-invariants into (what I call) “standard superfunctions”.

f(Z)×δ(C·η) -> {f(Z),C(Z)}

aka {bosonicPart,etaCoefficients}

Bosonic bits (trivial)

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

bosonicPart[exprn_]:=(exprn/.{R[x__]:>1/Times@@ab@@@Partition[{x},4,1,1]});

Fermionic parts (less trivial)

Note that each R-invariant gives it’s own row to C-matrix of η coefficients

rAmp[6,2]

R[cap[{2,3},{5,6,1}],3,4,5,cap[{5,6},{5,6,1}]]//.{R[x__]:>η/@{x}.ab@@@Partition[RotateLeft@{x},4,1,1],η[cap[{a_,b_},{cde__}]]:>η[a]ab[b,cde]+η[b]ab[cde,a]}
Normal[CoefficientArrays[%,η/@Range[6]][[-1]]]

What’s up with the caps?

η[cap[{2,3},{5,6,1}]]/.{η[cap[{a_,b_},{cde__}]]:>η[a]ab[b,cde]+η[b]ab[cde,a]}

Recall that the η’s are components of z’s, so caps are defined the same way:
(ab)⋂(cde):=z[a]⟨bcde⟩+z[b]⟨cdea⟩

etaCoefficients[n_][exprn_]:=
Block[{rList=Cases[exprn,R[__],{0,∞}],deltaArgs},
deltaArgs=rList//.{R[x__]:>η/@{x}.ab@@@Partition[RotateLeft@{x},4,1,1],η[cap[{a_,b_},{cde__}]]:>η[a]ab[b,cde]+η[b]ab[cde,a]};
Normal[CoefficientArrays[deltaArgs,η/@Range[n]][[-1]]]
]

Putting everything together

toSuperfunction[n_][exprn_]:={bosonicPart[exprn],etaCoefficients[n][exprn]};

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

Extracting component functions from superfunctions

Recall that each superfunction encodes ((Binomial[n,k])^4)(where 4 is 𝒩) bosonic components. Each of these are “bosonicPart×(four Det’s of etaCoeffs)”

We’re going to just build a single reference component extractor for the sake of checking things. This is related the so-call 
“split(-helicity (gluonic) amplitude)”

splitComponent[{bosonicPart_,etaCoefficients_}]:=
Block[{k=Length[etaCoefficients],minorList},
minorList=Range[k]&/@Range[4];
If[k==0,1,bosonicPart*Product[Det[etaCoefficients[[All,minor]]],{minor,minorList}]]
]

splitComponent/@superAmpRandom[8,1]

Evaluation and Verification

Strategy: two different methods of internal verification
1. numeric agreement among randomized formulae
2. check spurious pole cancellation (numerically)

Evaluation

We need to teach (very cautiously) Mathematica what an angle-bracket is. 
Recall that these are Det’s of momentum-twistor components. 

Note bene: because any unconstrained list of momentum-twistors encode
on-shell and momentum-conserving momenta, we may as well define a lot of them.

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

Teach mathematica what a 4-bracket is (temporarily!!)

Idea: don’t let evaluation persist outside of your needs

evaluate[exprn_]:=
Block[{output},
ab[x__Integer]:=ab[x]=Det[momentumTwistorComponents[[{x}]]];
ab[x___,cap[{a_,b_},{cde__}],y___]:=ab[x,a,y]ab[b,cde]+ab[x,b,y]ab[cde,a];
output=exprn;
ClearAll[ab];
output
]

evaluate[splitComponent/@superAmpRandom[10,3]];
Total[%]

Checking spurious pole cancellation

Recall that for planar amplitudes, the only physical poles are ones which set consecutive sets of momenta on-shell:

((((Subscript[p, a])+...+(Subscript[p, b-1])))^2)∝⟨a-1,a,b-1,b⟩->0

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

evaluate[physicalPoles[10]splitComponent/@superAmpRandom[10,2]];
IntegerQ/@%
Total[%%]//IntegerQ

Homework (expanded post-lecture)

Generalize the component extraction from the superfunctions

Recall that external states are specified by their “R-charges” (here, “R” refers to the (Subscript[SU(4), R])-symmetry).

(Subscript[(Subscript[(Subscript[(Subscript[(Subscript[|a⟩:=|a⟩, {}]^(+1))+|a⟩, {I}]^(+1/2))(Subscript[η, a]^I)+(Subscript[η, a]^I)(Subscript[η, a]^J)|a⟩, {I,J}]^0)+(Subscript[η, a]^I)(Subscript[η, a]^J)(Subscript[η, a]^K)|a⟩, {I,J,K}]^(-1/2))+(Subscript[η, a]^1)(Subscript[η, a]^2)(Subscript[η, a]^3)(Subscript[η, a]^4)|a⟩, {1,2,3,4}]^(-1))

Above, the R-charges of states are as follows: 
(+)-helicity gluon: {}
(+1/2)-helicity fermion indexed by I∈{1,...,4}: {I}
(0)-helicity scalar indexed by (the antisymmetric pair) I,J∈{1,...,4}: {I,J}
etc. 

So we can label the component function of a superfunction by the R-charges of the n external states. 
(Hint:) There must be k instances of each of the R-charge indices {1,...,4} among this list of data.

• Construct a function, rChargesToMinors[], which returns the list of 4(=𝒩) k×k minors (of the etaCoefficients matrix) which, when multiplied together (and when multiplied by the bosonicPart), give the component function involving that set of external states.

Nota bene: “momentum-twistor gluons” are not so-directly related to “momentum-space (aka actual) gluons”. This is because momentum-twistor external (“coherent”) wavefunctions are expanded in terms of η’s as opposed to (OverTilde[η])’s. (See below for more details.)

Uplift everything to momentum space (and extract real gluonic amplitudes)

From momentum-twistors to momenta (spinors)

Starting from momentum twistors, one constructs the data {λ,(OverTilde[λ])} as follows. 
(Subscript[z, a])=((Subscript[λ, a])|(Subscript[μ, a])). That is, the λ matrix is simply the top-two rows of the Z-matrix. 

From the μ’s, the (OverTilde[λ])’s are defined by 
(Subscript[(OverTilde[λ]), a]):=(Subscript[Q, ab])(Subscript[μ, b])
where 
(Subscript[Q, ab]):=(({{((⟨2n⟩)/(⟨n 1⟩⟨12⟩)), (1/(⟨12⟩)), 0, ⋯, ⋯, ⋯, (1/(⟨12⟩))}, {(1/(⟨12⟩)), ((⟨13⟩)/(⟨12⟩⟨23⟩)), (1/(⟨23⟩)), 0, ⋯, ⋯, 0}, {0, (1/(⟨23⟩)), ((⟨42⟩)/(⟨23⟩⟨34⟩)), (1/(⟨34⟩)), 0, ⋱, ⋮}, {⋮, 0, (1/(⟨34⟩)), ((⟨53⟩)/(⟨34⟩⟨45⟩)), (1/(⟨45⟩)), ⋱, ⋮}, {0, ⋮, ⋱, ⋱, ⋱, ⋱, 0}, {(1/(⟨n 1⟩)), 0, ⋯, ⋯, 0, (1/(⟨n-1 n⟩)), ((⟨1 n-1⟩)/(⟨n-1 n⟩⟨n 1⟩))}}))

(Recall that ⟨ab⟩:=Det[(Subscript[λ, a]),(Subscript[λ, b])] (which is in fact equal to  Det[momentumTwistorComponents[[{a,b},1;;2]]]).)

• Write a function which returns the above matrix for any n. 
• Check Q.λ=0—thus, (OverTilde[λ]) as defined this way is guaranteed to satisfy (OverTilde[λ]).λ=0.

From momenta (spinors) to momentum-twistors

From the (OverTilde[λ])’s, the μ’s may be defined by 
(Subscript[μ, a]):=(Subscript[(OverTilde[Q]), ab])(Subscript[(OverTilde[λ]), b])
where 
(Subscript[(OverTilde[Q]), ab]):=(({{0, 0, 0, ⋯, ⋯, ⋯, 0}, {0, 0, 0, 0, ⋯, ⋯, 0}, {0, ⟨12⟩, 0, 0, ⋱, ⋱, 0}, {0, ⟨24⟩, ⟨34⟩, 0, 0, ⋱, 0}, {0, ⟨25⟩, ⟨35⟩, ⟨45⟩, 0, ⋱, 0}, {⋮, ⋮, ⋮, ⋮, ⋱, ⋱, ⋮}, {0, ⟨2n⟩, ⟨3n⟩, ⟨4n⟩, ⋯, ⟨n-1 n⟩, 0}}))

• Write a function which returns the above matrix for any n.

Converting superfunctions

Make everything above more efficient!

For example, you can try to:
• Store values of rAmp permanently to memory (akin to termsInBCFW)
• Standardize 4-brackets so that all 4! different orderings need not be evaluated separately. (This improves evaluation speed.)
• Bake-in trivializations such as ⟨abc(cd)⋂(efg)⟩=⟨abcd⟩⟨efgc⟩.

Reference code (for the sake of comparison)

The following is how rAmp is defined in the package associated with the paper [1505.06137] (only slightly adapted for our present context).

bcfwBridgeRef[{{nL_,kL_},{nR_,kR_}}]:=Block[{n=nL+nR-2,a=nL,leftAmp,rightAmp,aHat,nHat},
aHat=cap[{a-1,a},{n-1,n,1}];
nHat=cap[{n-1,n},{a-1,a,1}];
leftAmp=rAmpRef[nL,kL]/.{R[x__]:>R@@({x}/.nL->aHat)};
rightAmp=(rAmpRef[nR,kR]/.{R[x__]:>R@@({x}/.Thread[Rule[Range[nR],Join[{aHat},Range[a,n-1],{nHat}]]])});
Times@@@Tuples[{leftAmp,{R[1,a-1,a,n-1,n]},rightAmp}]
];

ClearAll[rAmpRef];
rAmpRef[n_,k_]:=rAmpRef[n,k]=Which[k==0,{1},k>n-4,{},True,Join[rAmpRef[n-1,k],Join@@(bcfwBridgeRef/@bcfwPartitions[n,k])]];

Let’s compare speed (of (just) recursion) (just for fun :) )

Length[rAmpRef[12,4]]//Timing
Length[rAmp[12,4]]//Timing

Nota bene: the above speed comparison is really just about how quickly a formula for the n-point amplitude is found as a function of R-invariants. 

This is often not an important difference: once the fomula is obtained, it can be stored in memory for evaluation purposes (and the two algorithms rAmp and rAmpRef give the same expressions).