# Replace with non-commutative product

I'm trying to evaluate some scattering amplitudes terms in Mathematica using the spinor helicity formalism. For that, I need to turn products like $$\epsilon_i\cdot k_j$$ into spinor product which can be evaluated numerically. In order to do that, I'm using simple replacement rules in a function called subsp that has worked so far. However, now I'm encountering the following terms: $$\epsilon_i\cdot f_j\cdot f_k\cdot ...\cdot f_l \cdot \epsilon_p,$$

where the $$f$$ are defined as an outer product: $$f_i^{\mu\nu}=k_i^\mu\epsilon_i^\nu-\epsilon_i^\mu k_i^\nu$$. Regarding the code, my program firstly returns the terms like this: epsilon[i,poli].f[j,polj].f[k,polk].epsilon[l,poll] for example, where the different pol just indicate the polarization for each particle. Now, if I had products like $$\epsilon_i\cdot k_j$$, I would just make replacements like

epsilon[a_, P].k[b_Integer] -> epsPk[a, b, q_a],

where epsPk[a, b, q_a] is a function that knows how to evaluate such product with spinors (the P just indicates positive polarization). For the $$f$$, I tried doing

f[i_, poli_] :> k[i]epsilon[i,poli]-epsilon[i,poli]k[i],

but obviously this evaluates to zero. I also tried replacing products like $$\epsilon_i\cdot f_j$$ with

epsilon[i_, poli_].f[j_, polj_] -> epsilon[i, poli].k[j]*epsilon[j, polj] - epsilon[i, poli].epsilon[j, polj]*k[j],

but this also doesn't work. The problem is that longer chains get blocked by the priority of Times over Dot (I think that's what the problem is). For example, if I type:

subsp[epsilon[3, p].f[2, p].k[4]]

I get the following output:

-(k[2]spb[2, 3]).k[4]/spa[3, 2] + Sqrt[2] (epsilon[2, p]spa[4, 2]spb[2, 3]).k[4]/spa[4, 3]

As you can see, the function has succesfully contracted epsilon[3,p] with the fisrt half of f[2,p] (notice the spinor terms like spa[3, 2] or spb[2, 3]), but fails to do so with the right half of f[2,p] and k[4] because it is blocked by the parenthesis (more precisely, I would need it to further contract k[2].k[4] in the negative term and epsilon[2,p].k[2] in the positive one).

In conclusion, I need to be able to successfully define either a non-vanishing f[i,poli] or an arbitrarily long chain of products epsilon[i,poli].f[j,polj].....f[k,polk].epsilon[l,poll] with the replacement rules. If anyone could help me with this issue, I would reall appreciate it. Thanks in advance.

Edit: I've been requested to provide more insight into what my code looks like and what I'm aiming. What I'm trying to calculate specifically are something called BCJ numerators, which I obtain from a set of tree diagrams with another function that isn't important to this issue. One possible output that can serve as the starting point to illustrate the problem is:

epsilon[3,p].epsilon[5,m] k[3].k[{1,2}]epsilon[1,p].f[3,p].epsilon[5,m]

The only new thing apart from what I've already explained are lists being arguments of k. This represents a sum of momenta, e.g. k[{1,2}]=k[1]+k[2] (which I've succesfully implemented in the function below). I would like to keep this output format intact. Now, I want subsp to convert this into a product of spinor functions spa and spb. This is howI do it:

subsp[expr_, n_] := Module[{a, b, c, pol, pola, polb}, rules = {s[a_, b_] -> spa[a, b]*spb[b, a], epsPP[a_, b_, qa_, qb_] -> spa[qa, qb]*spb[b, a]/(spa[qa, a]*spa[qb, b]), epsPM[a_, b_, qa_, qb_] -> spa[qa, b]*spb[qb, a]/(spa[qa, a]*spb[qb, b]), epsMP[a_, b_, qa_, qb_] -> spb[qa, b]*spa[qb, a]/(spb[qa, a]*spa[qb, b]), epsMM[a_, b_, qa_, qb_] -> spb[qa, qb]*spa[b, a]/(spb[qa, a]*spb[qb, b]), epsPk[a_, b_, qa_] -> Sqrt[2]*spa[qa, b]*spb[b, a]/spa[qa, a], epsMk[a_, b_, qa_] -> Sqrt[2]*spb[qa, b]*spa[b, a]/spb[qa, a], epsilon[b_, pol_].k[a : _List] :> TensorExpand[epsilon[b, pol].Total[ Map[k, a]]], epsilon[a_, p].epsilon[b_, p] -> epsPP[a, b, Mod[a + 1, n, 1], Mod[b + 1, n, 1]], epsilon[a_, p].epsilon[b_, m] -> epsPM[a, b, Mod[a + 1, n, 1], Mod[b + 1, n, 1]], epsilon[a_, m].epsilon[b_, p] -> epsMP[a, b, Mod[a + 1, n, 1], Mod[b + 1, n, 1]], epsilon[a_, m].epsilon[b_, m] -> epsMM[a, b, Mod[a + 1, n, 1] + 1, Mod[b + 1, n, 1]], epsilon[a_, p].k[b_Integer] -> epsPk[a, b, Mod[a + 1, n, 1]], epsilon[a_, m].k[b_Integer] -> epsMk[a, b, Mod[a + 1, n, 1]], k[a_Integer]. k[b_List] :> TensorExpand[k[a].Total[Map[k, b]]], k[a_Integer].k[b_Integer] -> 1/2*s[a, b], f[b_, polb_] :> k[b] ** epsilon[b, polb] - epsilon[b, polb] ** k[b], expf[a_] -> 1 + a};expr2 = FixedPoint[Expand[TensorExpand[ReplaceRepeated[#, rules]]] &, expr]];

I know this has some inefficient definitions like the redundant intermediate functions epsPP or epskM, but they are there because they are useful in some other applications. Using this function, the DESIRED output of the previous example would be:

(Sqrt[2] spa[3, 5] spa[4, 5] spb[6, 3]^2)/(spa[4, 3] spb[6, 5]^2) - (Sqrt[2] spa[2, 4] spa[3, 5] spb[3, 1] spb[6, 3])/(spa[2, 1] spa[4, 3] spb[6, 5]),


but the problem I mentioned with the f functions instead gives me:

-((epsilon[1, p].epsilon[3, p] ** k[3].epsilon[5, m] spa[1, 3] spa[4,5] spb[3, 1] spb[6, 3])/(2 spa[4, 3] spb[6, 5])) + (epsilon[1, p].k[3] ** epsilon[3, p].epsilon[5, m] spa[1, 3] spa[4,5] spb[3, 1] spb[6, 3])/(2 spa[4, 3] spb[6, 5]) - (epsilon[1, p].epsilon[3, p] ** k[3].epsilon[5, m] spa[2, 3] spa[4, 5] spb[3, 2] spb[6, 3])/(2 spa[4, 3] spb[6, 5]) + (epsilon[1, p].k[3] ** epsilon[3, p].epsilon[5, m] spa[2, 3] spa[4, 5] spb[3, 2] spb[6, 3])/(2 spa[4, 3] spb[6, 5])


i.e. there are still epsilon and k functions to be contracted.

• I'd suggest NonCommutativeMultiply (**) for your products - as the name suggests, it's not commutative, and it also doesn't have any built-in meaning that would lead to unwanted evaluation Commented Sep 6, 2022 at 13:45
• I tried that. However, the problem persists: NonCommutativeMultiply has priority with respect to Dot so, even after TensorExpand, the output looks like -epsilon[3, p].epsilon[2, p] ** k[2].k[4] + epsilon[3, p].k[2] ** epsilon[2, p].k[4], with no further evaluation Commented Sep 6, 2022 at 13:59
• Can you please provide a fully working sample of your current code with the current & the desired output? E.g. you mention subsp, but it's not clear how it is defined Commented Sep 6, 2022 at 14:29
• @LukasLang I've added more details in the post and provided an example. Commented Sep 6, 2022 at 20:13

This is not yet a complete answer, but I hope it's a start: It manages to contract all the epsilon factors, but the result is not what you mention as expected output. Unfortunately, I don't understand your notation enough to tell where the problem lies.

subsp[expr_, n_] :=
Module[{rules},
rules = {s[a_, b_] :> spa[a, b]*spb[b, a],
epsPP[a_, b_, qa_, qb_] :>
spa[qa, qb]*spb[b, a]/(spa[qa, a]*spa[qb, b]),
epsPM[a_, b_, qa_, qb_] :>
spa[qa, b]*spb[qb, a]/(spa[qa, a]*spb[qb, b]),
epsMP[a_, b_, qa_, qb_] :>
spb[qa, b]*spa[qb, a]/(spb[qa, a]*spa[qb, b]),
epsMM[a_, b_, qa_, qb_] :>
spb[qa, qb]*spa[b, a]/(spb[qa, a]*spb[qb, b]),
epsPk[a_, b_, qa_] :> Sqrt[2]*spa[qa, b]*spb[b, a]/spa[qa, a],
epsMk[a_, b_, qa_] :> Sqrt[2]*spb[qa, b]*spa[b, a]/spb[qa, a],
epsilon[b_, pol_] ** Dot ** k[a : _List] :>
epsilon[b, pol] ** Dot ** Total[Map[k, a]],
epsilon[a_, p] ** Dot ** epsilon[b_, p] :>
epsPP[a, b, Mod[a + 1, n, 1], Mod[b + 1, n, 1]],
epsilon[a_, p] ** Dot ** epsilon[b_, m] :>
epsPM[a, b, Mod[a + 1, n, 1], Mod[b + 1, n, 1]],
epsilon[a_, m] ** Dot ** epsilon[b_, p] :>
epsMP[a, b, Mod[a + 1, n, 1], Mod[b + 1, n, 1]],
epsilon[a_, m] ** Dot ** epsilon[b_, m] :>
epsMM[a, b, Mod[a + 1, n, 1] + 1, Mod[b + 1, n, 1]],
epsilon[a_, p] ** Dot ** k[b_Integer] :>
epsPk[a, b, Mod[a + 1, n, 1]],
epsilon[a_, m] ** Dot ** k[b_Integer] :>
epsMk[a, b, Mod[a + 1, n, 1]],
k[b_Integer] ** Dot ** epsilon[a_, p] :>
epsPk[a, b, Mod[a + 1, n, 1]],
k[b_Integer] ** Dot ** epsilon[a_, m] :>
epsMk[a, b, Mod[a + 1, n, 1]],
k[a_Integer] ** Dot ** k[b_List] :>
k[a] ** Dot ** Total[Map[k, b]],
k[a_Integer] ** Dot ** k[b_Integer] :> 1/2*s[a, b],
f[b_, polb_] :>
k[b] ** epsilon[b, polb] - epsilon[b, polb] ** k[b],
expf[a_] :> Identity + a,
Identity ** Dot ** a_ | a_  ** Dot ** Identity :> a,
(m_?NumericQ a_) ** b_ :> a ** b,
a_ ** (m_?NumericQ b_) :> m a ** b,
a_ ** (b_ / c_) :> a ** b / c,
(a_ / c_) ** b_ :> a ** b / c};
FixedPoint[
ReplaceRepeated[# /. m_NonCommutativeMultiply :> Distribute[m],
rules] &, expr] /. {NonCommutativeMultiply -> Times}];

subsp[
epsilon[3, p] ** Dot ** epsilon[5, m] ** k[3] ** Dot ** k[{1, 2}] **
epsilon[1, p] ** Dot ** f[3, p] ** Dot ** epsilon[5, m], 100] // FullSimplify
(* (spa[4,
5] (-spa[2, 4] spa[3, 5] +
spa[2, 3] spa[4, 5]) (spa[3, 1] spb[1, 3] +
spa[3, 2] spb[2, 3]) spb[3, 1] spb[6, 3]^2)/(Sqrt[2]
spa[2, 1] spa[4, 3]^2 spb[6, 5]^2) *)

subsp[
epsilon[1, p] ** Dot ** f[3, p] ** Dot ** epsilon[5, m], 100] // FullSimplify // Expand
(* -((Sqrt[2] spa[2, 4] spa[3, 5] spb[3, 1] spb[6, 3])/(
spa[2, 1] spa[4, 3] spb[6, 5])) + (
Sqrt[2] spa[2, 3] spa[4, 5] spb[3, 1] spb[6, 3])/(
spa[2, 1] spa[4, 3] spb[6, 5]) *)


• Made proper use of Module to localize rules instead of the (already localized) pattern names
• Replaced everything with RuleDelayed to help the syntax highlighter
• Replaced all relevant products with **. The products that indicate dot products are written as a**Dot**b. This prevents any issues with incorrect grouping between e.g. . and ** products that you have seen (thanks to the Flat attribute of **) while still allowing us to distinguish "normal" products from dot products.
• I added rules for k**Dot**epsilon style dot products (these might be wrong, currently they are the same as the corresponding epsilon**Dot**k products)
• I removed all the useless Expand and TensorExpand calls
• I added two rules (m_?NumericQ a_) ** b_ :> m a ** b, a_ ** (m_?NumericQ b_) :> m a ** b to factor out numeric prefactors that interfere with the other rules
• Similarly, I added two rules a_ ** (b_ / c_) :> a ** b / c, (a_ / c_) ** b_ :> a ** b / c that move denominators to the end of the prosuct, such that they are out of the way of the other rules
• I added a call to Distribute to distribute the ** over the +.
• I replace ** with the normal * in the last step to allow the final result to properly simplify.
• I replaced the 1 in expf with Identity to indicate the identity tensor, and added the relevant contration rules

I hope this helps you in some way, even though the result is not yet what you expect (as mentioned).

With all that being said, I am not convinced this is the best route forward: It seems to me that a formalism capturing the indices being contracted (the $$\mu$$ and $$\nu$$ in your example) would be way easier. That removes the issues with products that should not commute and automatically gives you the distributive law and factoring out of numeric prefactors etc.

• After some minor tweeking, this definitely works better than what I had, but there are still some problematic cases. For example, if I type this: subsp[epsilon[2, p]**Dot**k[{1}]**epsilon[3, m]**Dot**k[{1, 2}],4], I get this result: Sqrt[2]((spa[3, 1] spb[1, 2])/spa[3, 2])**epsilon[3, m]**Dot**(k[1] + k[2]), which is still not completely reduced. It seems like the Distribute isn't working as it should. Commented Sep 7, 2022 at 10:53
• @Marcosko I have edited the answer with a fix for the Distribute issue. The issue was that the expression was something like a*(b**c), meaning the Distribute didn't see the **. The new version simply applies Distribute manually. As for the other tweaks: In the interest of making the answer fully correct & working, what changes were required? (Please feel free to suggest an edit to the answer if that's easier than describing the fixes in the comments) Commented Sep 7, 2022 at 11:35
• Don't worry too much about the tweeks, they were just to adapt the function to additional stuff I didn't include in the post in order to not complicate things (I probably should have phrased my comment differently). Anyways, there are two remaining problems now: first, the function fails when there are fractions, sine it can't differentiate numerators and denominators faithfully when ussing ** and **Dot**. For example, trying to calculate subsp[epsilon[1, p]**Dot**(f[2, p]/epsilon[2, p]**Dot**k[{1}])**Dot**(f[3, m]/epsilon[3, m]**Dot**k[{1, 2}])**Dot** epsilon[4, m], 7] fails. Commented Sep 7, 2022 at 12:24
• @Marcosko what is the intended meaning of the division? The same as for multiplication, but with a 1/epsilon factor? Commented Sep 7, 2022 at 12:27
• The second problem comes from the presence of the exponential expf which I expand as expf[a_] :> 1 + a. This leads to some terms having things like **Dot**1**Dot** with which the function doesn't know how to deal. Have you any idea of how it could be fixed? Commented Sep 7, 2022 at 12:27

One could define

dots[x__,f[b_,polb_],y__]:=( dots[x,k[b]]*dots[epsilon[b,polb],y]
-dots[x,epsilon[b,polb]]*dots[k[b],y])//Expand;


For example,

dots[epsilon[1,p],f[2,p],f[3,m],k[4]]


returns

dots[epsilon[1,p],k[2]] dots[epsilon[2,p],k[3]] dots[epsilon[3,m],k[4]]
-dots[epsilon[1,p],epsilon[2,p]] dots[epsilon[3,m],k[4]] dots[k[2],k[3]]
-dots[epsilon[1,p],k[2]] dots[epsilon[2,p],epsilon[3,m]] dots[k[3],k[4]]
+dots[epsilon[1,p],epsilon[2,p]] dots[k[2],epsilon[3,m]] dots[k[3],k[4]]


which is a sum of products of simple terms of the form dots[k[_],epsilon[_,_]] and similar, which can then be rewritten in terms of spa and spb.

• Yes, but each element x, f, y are not connected in the original output (and it would be very non-trivial to do so), i.e. they can't appear simply as arguments of the dots` function. Commented Sep 8, 2022 at 5:20