Bug with FCFADiracChainJoin

[Turgon-Aran-Gondolin]

• Your Mathematica version

13.0.1 for Microsoft Windows (64-bit) (January 28, 2022)

• Your FeynCalc version

9.3.1

• Did you try to reinstall FeynCalc (stable version) using the automatic installer to make sure that you have the latest bugfixes?

No, but I did check with the latest master branch.

• Does your Mathematica initialization file contain statements that might influence the behavior of FeynCalc? Sometimes external packages may modify `init.m` in unusual ways, causing troubles for other codes.

Nothing.

• Please provide a minimal working example that illustrates the problem and works on a fresh kernel. The example should be provided either by writing the code (as `InputForm`!) directly in the issue text or by attaching a Mathematica notebook. Please do not post code samples as screenshots, PDF files etc.: Those essentially require us to retype everything by hand, which is annoying and also time consuming. Please explain the difference between the current behavior and the expected behavior.

DiracChain[Spinor[-Momentum[p1, D], mqu, 1], 
  DiracIndex[Dir1]] DiracChain[Spinor[-Momentum[p2, D], mqu, 1], 
  DiracIndex[Dir2]] DiracChain[DiracGamma[LorentzIndex[Lor2, D], D],
   DiracIndex[Dir5], DiracIndex[Dir1]] DiracChain[
  DiracGamma[LorentzIndex[Lor3, D], D], DiracIndex[Dir2], 
  DiracIndex[Dir7]] DiracChain[
  mqu + DiracGamma[Momentum[-l - p2, D], D], DiracIndex[Dir7], 
  DiracIndex[Dir8]] DiracChain[
  mqu + DiracGamma[Momentum[-k - k1 - k2 + l + p2, D], D], 
  DiracIndex[Dir5], DiracIndex[Dir6]]  DiracChain[
  DiracGamma[LorentzIndex[LCdummy2, D], D] . DiracGamma[6], 
  DiracIndex[Dir8], DiracIndex[Dir6]]] // FCE

output:
-Spinor[-Momentum[p1, D], mqu, 1] . (-GAD[Lor2]) . (mqu + GSD[-k - k1 - k2 + l + p2]) . GA[6] . GAD[LCdummy2] . (mqu - GSD[-l - p2]) . GAD[Lor3] . Spinor[Momentum[p2, D], mqu, 1]

The Dirac chain should be in reversed order.

[vsht]

The task of FCFADiracChainJoin is merely to create a proper chain from the output of FeynArts.
The ordering of the spinors can be pretty much arbitrary at this stage.

However, you can always reverse the given chain using SpinorChainTranspose. Using suitable patterns via the Select option you can explicitly specify, which chains must be reversed and, which not. For example

-Spinor[-Momentum[p1, D], mqu, 
    1] . (-GAD[Lor2]) . (mqu + GSD[-k - k1 - k2 + l + p2]) . GA[6] . 
   GAD[LCdummy2] . (mqu - GSD[-l - p2]) . GAD[Lor3] . 
   Spinor[Momentum[p2, D], mqu, 1] // 
 SpinorChainTranspose[#, FCE -> True] &

returns

Spinor[-Momentum[p2, D], mqu, 1] . GAD[Lor3] . (mqu + GSD[-l - p2]) . 
 GAD[LCdummy2] . GA[6] . (mqu - GSD[-k - k1 - k2 + l + p2]) . 
 GAD[Lor2] . Spinor[Momentum[p1, D], mqu, 1]

Does this solve the problem?

[Turgon-Aran-Gondolin]

Thanks, SpinorChainTranspose provides a temporary solution. However, since FCFADiracChainJoin is called in the last step of FCFAConvert (which is where I initially encountered this problem), it would be better if the order is properly dealt with by FCFADiracChainJoin. Calling SpinorChainTranspose is especially tedious if there're multiple spinor chains in the amplitude, and one must look at each individual diagram to determine if this extra step is needed.

I'd argue that the ordering of the spinors is not arbitrary, in the sense that the indices of DiracChain should not be orderless. In the example I gave, the first and the last Dirac matrice read DiracChain[ DiracGamma[LorentzIndex[Lor3, D], D], DiracIndex[Dir2], DiracIndex[Dir7]] and DiracChain[DiracGamma[LorentzIndex[Lor2, D], D], DiracIndex[Dir5], DiracIndex[Dir1]]. Assuming the Dirac indices are not orderless, the first index should be Dir2, and the last Dir1. So the bilinear should start with a spinor with index Dir2, and end with Dir1.

[vsht]

I'd argue that the ordering of the spinors is not arbitrary, in the sense that the indices of DiracChain should not be orderless. In the example I gave, the first and the last Dirac matrice read DiracChain[ DiracGamma[LorentzIndex[Lor3, D], D], DiracIndex[Dir2], DiracIndex[Dir7]] and DiracChain[DiracGamma[LorentzIndex[Lor2, D], D], DiracIndex[Dir5], DiracIndex[Dir1]]. Assuming the Dirac indices are not orderless, the first index should be Dir2, and the last Dir1. So the bilinear should start with a spinor with index Dir2, and end with Dir1.

The indices are not orderless, but the choice of the matrix that will be complex conjugated first is indeed pretty much arbitrary.
You would like it to start with {Dir2,Dir7} {Dir2,Dir5}, but starting with say {Dir5,Dir1} {Dir5,Dir6} or {Dir8,Dir6} {Dir5,Dir6} is equally legitimate. Ultimately, it is the way how the rules are internally sorted by Mathematica's evaluator which decides which pair of matrices will be treated first.

Forcing it to start with Dir2 would require to build all possible matrix pairs first and identify the lexicographically
lowest index. This would significantly complicate the algorithm without really gaining much.

Calling SpinorChainTranspose is especially tedious if there're multiple spinor chains in the amplitude, and one must look at each individual diagram to determine if this extra step is needed.

Why? You can just define the ordering you want via the Select option and then apply it to all diagrams. If the ordering is already there, nothing would happen. Otherwise, the spinors get reordered accordingly. It is just a single step needed at the beginning of the calculation (even before doing any evaluations/simplifications).

I use a similar approach e,g. to fix the absolute signs of the amplitudes to my liking, since those are a pure convention. Just like the ordering of the spinors in this case.

Apart from that, now I remember that there is actually an option First in FCFADiracChainJoin for choosing the spinor that
should come first in the final chain. It still uses SpinorChainTranspose behind the scenes, but perhaps it suits your needs in a better way. The option was, in fact, broken, since I never used it before, but this should be now fixed.

So if you use FCFAConvert with FCFADiracChainJoin set to False, and then run FCFADiracChainJoin manually
with the option First, this should do the job

chain = DiracChain[Spinor[-Momentum[p1, D], mqu, 1], 
   DiracIndex[Dir1]] DiracChain[Spinor[-Momentum[p2, D], mqu, 1], 
   DiracIndex[Dir2]] DiracChain[DiracGamma[LorentzIndex[Lor2, D], D], 
   DiracIndex[Dir5], DiracIndex[Dir1]] DiracChain[
   DiracGamma[LorentzIndex[Lor3, D], D], DiracIndex[Dir2], 
   DiracIndex[Dir7]] DiracChain[
   mqu + DiracGamma[Momentum[-l - p2, D], D], DiracIndex[Dir7], 
   DiracIndex[Dir8]] DiracChain[
   mqu + DiracGamma[Momentum[-k - k1 - k2 + l + p2, D], D], 
   DiracIndex[Dir5], DiracIndex[Dir6]] DiracChain[
   DiracGamma[LorentzIndex[LCdummy2, D], D] . DiracGamma[6], 
   DiracIndex[Dir8], DiracIndex[Dir6]];

FCFADiracChainJoin[chain, First -> {Spinor[-Momentum[p2, D], mqu, 1]}]
(*or*)
FCFADiracChainJoin[chain, First -> {Spinor[-Momentum[p1, D], mqu, 1]}]

Does this help?

[Turgon-Aran-Gondolin]

What I'm proposing is actually this: first find out what's the indices of the spinors, say {Dir1, Dir2}; then search for the matrice with these two indices, say the result is {{Dir5, Dir1}, {Dir2, Dir7}}. Then because Dir2 is the first index in its own pair, Dir1 is the last, we can say Dir2 is related to the first spinor, and Dir1 the last.

Thanks for sharing with me the First option, I didn't know that. And you are absolutely right: getting the order sorted just requires some simple additional steps. I just thought it would be great if those are sorted automatically.

[vsht]

What you suggest is a prescription that is, in principle, as good as any other, but technically, it is more expensive
performance wise. The way FCFADiracChainJoin works is to assemble the matrix chain first and then
connect the spinors in the very last step, simultaneously fixing the relative fermion sign.

Starting with the spinors and using open indices to determine, which matrices should be CCed
first, adds a layer of complexity on top of that, where the benefits are, in my view, rather opinion-based:
Next time someone else would complain that your prescription doesn't suit him/her and suggest
another ordering method.

I think that with the First option and explicit SpinorChainTranspose there is already enough flexibility for those
who care about the ordering. Those who just need to square the amplitude wouldn't care anyway.

Having said that, you can of course write and use your own version of FCFADiracChainJoin that does exactly
what you want. The source code is rather short and the algorithm should be transparent enough.

[Turgon-Aran-Gondolin]

I tried the First option, some of the amplitudes are fine, but this failed:

FCFADiracChainJoin[
 DiracChain[Spinor[Momentum[k1, D], mle, 1], 
   DiracIndex[Dir3]] DiracChain[Spinor[Momentum[k2, D], mle, 1], 
   DiracIndex[Dir4]] DiracChain[Spinor[-Momentum[p1, D], mqu, 1], 
   DiracIndex[Dir1]] DiracChain[Spinor[-Momentum[p2, D], mqu, 1], 
   DiracIndex[Dir2]] DiracChain[DiracGamma[LorentzIndex[Lor2, D], D], 
   DiracIndex[Dir5], DiracIndex[Dir1]] DiracChain[
   mqu + DiracGamma[Momentum[-k - k1 - k2 + l + p2, D], D], 
   DiracIndex[Dir5], DiracIndex[Dir6]] DiracChain[
   DiracGamma[LorentzIndex[LCdummy1, D], D] . DiracGamma[6], 
   DiracIndex[Dir3], DiracIndex[Dir4]] DiracChain[
   DiracGamma[LorentzIndex[LCdummy2, D], D] . DiracGamma[6], 
   DiracIndex[Dir2], DiracIndex[Dir6]], 
 First -> {Spinor[-Momentum[p2, D], mqu, 1], 
   Spinor[Momentum[k1, D], mle, 1]}]

Do you know what went wrong?

I still don't understand why there could be multiple descriptions of ordering. I believe the ordering is fixed once the matrice with ordered indices are given. After contracting the matrice, the chain of the matrice is left with two open indices, and their positions are fixed. Then the way of adding spinors to the chain is also fixed. Can you provide an example of alternative descriptions?

[vsht]

Do you know what went wrong?

This should be also now fixed, thanks.

I still don't understand why there could be multiple descriptions of ordering. I believe the ordering is fixed once the matrice with ordered indices are given. After contracting the matrice, the chain of the matrice is left with two open indices, and their positions are fixed. Then the way of adding spinors to the chain is also fixed. Can you provide an example of alternative descriptions?

If you ignore the spinors for a moment and start contracting matrices, you will encounter multiple products of matrices that have their indices written in a "wrong" way so that you can't contract them. Using charge conjugation you can rewrite them as follows

$A_{ai} B_{bi} \to (A. C B^T C^-1)_{ab} $

$A_{ia} B_{ib} \to (C A^T C^-1)_{ab}$

However, you are completely free to choose the $A$ and $B$ matrices you start with. Different choices may lead to the final chain coming up as $X_{ab}$ or $X_{ba}$.

Then you attach the external spinors $\phi_a$ and
$\chi_{b}$ and depending on whether you got $X_{ab}$ or $X_{ba}$, a different ordering arises.

[Turgon-Aran-Gondolin]

Thanks, you gave some very compelling examples. However, the fermion loop is the only scenario I'm aware of with this kind of "flipped" spinor indices. Are there some new physics models that have this kind of structure?

Another issue is the fermion sign:

FCFADiracChainJoin[
 DiracChain[Spinor[Momentum[k1, D], mle, 1], 
   DiracIndex[Dir3]] DiracChain[Spinor[Momentum[k2, D], mle, 1], 
   DiracIndex[Dir4]] DiracChain[Spinor[-Momentum[p1, D], mqu, 1], 
   DiracIndex[Dir1]] DiracChain[Spinor[-Momentum[p2, D], mqu, 1], 
   DiracIndex[Dir2]] DiracChain[DiracGamma[LorentzIndex[Lor2, D], D], 
   DiracIndex[Dir5], DiracIndex[Dir1]] DiracChain[
   mqu + DiracGamma[Momentum[-k - k1 - k2 + l + p2, D], D], 
   DiracIndex[Dir5], DiracIndex[Dir6]] DiracChain[
   DiracGamma[LorentzIndex[LCdummy1, D], D] . DiracGamma[6], 
   DiracIndex[Dir3], DiracIndex[Dir4]] DiracChain[
   DiracGamma[LorentzIndex[LCdummy2, D], D] . DiracGamma[6], 
   DiracIndex[Dir2], DiracIndex[Dir6]], 
 First -> {Spinor[-Momentum[p2, D], mqu, 1], 
   Spinor[Momentum[k1, D], mle, 1]}]

and

FCFADiracChainJoin[
  DiracChain[Spinor[Momentum[k1, D], mle, 1], 
    DiracIndex[Dir3]] DiracChain[Spinor[Momentum[k2, D], mle, 1], 
    DiracIndex[Dir4]] DiracChain[Spinor[-Momentum[p1, D], mqu, 1], 
    DiracIndex[Dir1]] DiracChain[Spinor[-Momentum[p2, D], mqu, 1], 
    DiracIndex[Dir2]] DiracChain[DiracGamma[LorentzIndex[Lor2, D], D],
     DiracIndex[Dir5], DiracIndex[Dir1]] DiracChain[
    mqu + DiracGamma[Momentum[-k - k1 - k2 + l + p2, D], D], 
    DiracIndex[Dir5], DiracIndex[Dir6]] DiracChain[
    DiracGamma[LorentzIndex[LCdummy1, D], D] . DiracGamma[6], 
    DiracIndex[Dir3], DiracIndex[Dir4]] DiracChain[
    DiracGamma[LorentzIndex[LCdummy2, D], D] . DiracGamma[6], 
    DiracIndex[Dir2], DiracIndex[Dir6]]] /. 
 Spinor[-Momentum[p1, D], mqu, 1] . a__ . 
   Spinor[Momentum[p2, D], mqu, 1] :> 
  SpinorChainTranspose[
   Spinor[-Momentum[p1, D], mqu, 1] . a . 
    Spinor[Momentum[p2, D], mqu, 1]]

gives different signs. By checking the gauge invariance (in my own problem with other diagrams), the latter seems to have the correct sign.

[vsht]

The second code block seems to have a syntax error and doesn't produce a closed chain, or am I missing something?

[Turgon-Aran-Gondolin]

My apologies. There's a FCFADiracChainJoin[ missing in the second code block.

[vsht]

Thanks. This is indeed a very stupid bug from my side. The momentum-dependent ordering should
of course remain the same, since the transposition is done via a separate function.

I hope, this should be now fixed in both versions.

[Turgon-Aran-Gondolin]

Thanks. I'm closing this issue, as I'm not seeing any problems at the moment.

[Turgon-Aran-Gondolin]

I encountered another problem with the First option.

Consider the following DiracChain object:

dcobj = DiracChain[Spinor[Momentum[k1, D], MLE, 1], 
   DiracIndex[Dir3]] DiracChain[Spinor[Momentum[k2, D], MLE, 1], 
   DiracIndex[Dir4]] DiracChain[Spinor[-Momentum[p1, D], MQU, 1], 
   DiracIndex[Dir1]] DiracChain[Spinor[-Momentum[p2, D], MQU, 1], 
   DiracIndex[Dir2]] DiracChain[DiracGamma[LorentzIndex[Lor1, D], D], 
   DiracIndex[Dir2], DiracIndex[Dir5]] DiracChain[
   DiracGamma[LorentzIndex[Lor2, D], D], DiracIndex[Dir7], 
   DiracIndex[Dir4]] DiracChain[
   MLE + DiracGamma[Momentum[k2 - l, D], D], DiracIndex[Dir7], 
   DiracIndex[Dir8]] DiracChain[
   MQU + DiracGamma[Momentum[l - p2, D], D], DiracIndex[Dir5], 
   DiracIndex[Dir6]] DiracChain[
   DiracGamma[LorentzIndex[Ind37202, D], D] . DiracGamma[6], 
   DiracIndex[Dir3], DiracIndex[Dir8]] DiracChain[
   DiracGamma[LorentzIndex[Ind37202, D], D] . DiracGamma[6], 
   DiracIndex[Dir6], DiracIndex[Dir1]]

Doing FCFADiracChainJoin with a First option gives Evaluation of isolated objects failed error:

FCFADiracChainJoin[dcobj, First -> {Spinor[-Momentum[k1, D], MLE, 1]}]

I tested it with your latest hotfix branch. The expected behavior is to have a ubar(k1) v(k2) bilinear, but the default without First produces ubar(k2) v(k1).

Also, the default Select option of SpinorChainTranspose cannot handle ubar v type bilinear. Of course I can add it manually (presumably with a {SpinorUBarD[_, _], SpinorVD[_, _]}?), but I just want to ask if there's any reason you didn't add it in the first place?

[vsht]

I tested it with your latest hotfix branch. The expected behavior is to have a ubar(k1) v(k2) bilinear, but the default without First produces ubar(k2) v(k1).

Thanks, this should be now fixed in master.

Also, the default Select option of SpinorChainTranspose cannot handle ubar v type bilinear. Of course I can add it manually (presumably with a {SpinorUBarD[_, ], SpinorVD[, _]}?), but I just want to ask if there's any reason you didn't add it in the first place?

Because when this function is applied standalone (not inside FermionSpinSum) you usually want to get the same ordering in all spinor chains, like with

a SpinorVBar[p1, m1] . GS[p] . SpinorU[p2, m2] +  b SpinorVBar[p1, m1] . GS[p] . SpinorU[p2, m2]
% // SpinorChainTranspose

It makes not so much sense to just transpose everything. The default behavior is also documented and using the Select option
everyone can choose whatever they like.

[Ichthyocentaur]

• Mathematica version
  13.1 on Windows
• FeynCalc version
  10.0.0
• Did you try to reinstall FeynCalc (stable version) using the automatic installer to make sure that you have the latest bugfixes?
  Yes
• Does your Mathematica initialization file contain statements that might influence the behavior of FeynCalc? Sometimes external packages may modify init.m in unusual ways, causing troubles for other codes.
  No
• MWE
  I have attached my FeynRules model file and Mathematica notebook, and when I use the command FCFAConvert, I get the error

"Error! FCFADiracChainJoin has encountered a fatal problem and must \
abort the computation. The problem reads: \
\!\(TraditionalForm\`\"Evaluation of isolated objects failed.\"\)"

Model File

M$ModelName = "FermionEFT";

M$Parameters = {
    m == {ParameterType -> Internal},
    c1 == {ParameterType -> External},
    c2 == {ParameterType -> External},
    c3 == {ParameterType -> External}
};

M$ClassesDescription = {
    F[1] == {
        ClassName -> psi,
        Mass -> m,
        PropagatorType -> Straight,
        PropagatorArrow -> Forward,
        SelfConjugate -> False
    }
};

L4 = I psibar.Ga[mu].del[psi, mu] - m psibar.psi;

L6 = c1 Module[ {i,j}, psibar[i].psi[i].psibar[j].psi[j] ] + c2 Module[ {i,j,k,l}, Ga[mu, i, j] Ga[mu, k, l] psibar[i].psi[j].psibar[k].psi[l] ];

L8 = c3 Module[ {i,j}, psibar[i].del[psi[i], mu].psibar[j].del[psi[j], mu] ];

LFermionEFT = L4 + L6;

Mathematica file

description="ee -> ee, Fermion EFT, tree";
If[ $FrontEnd === Null,
	$FeynCalcStartupMessages = False;
	Print[description];
];
If[ $Notebooks === False,
	$FeynCalcStartupMessages = False
];
$LoadAddOns={"FeynArts","FeynHelpers"};
<<FeynCalc`
$FAVerbose = 0;

FCCheckVersion[9,3,1];
MakeBoxes[p1,TraditionalForm]:="\!\(\*SubscriptBox[\(p\), \(1\)]\)";
MakeBoxes[p2,TraditionalForm]:="\!\(\*SubscriptBox[\(p\), \(2\)]\)";
MakeBoxes[p3,TraditionalForm]:="\!\(\*SubscriptBox[\(p\), \(3\)]\)";
MakeBoxes[p4,TraditionalForm]:="\!\(\*SubscriptBox[\(p\), \(4\)]\)";
MakeBoxes[c1,TraditionalForm]:="\!\(\*SubscriptBox[\(c\), \(1\)]\)";
FAPatch[PatchModelsOnly->True];

diags = InsertFields[CreateTopologies[0, 2 -> 2,
		ExcludeTopologies->{}],
		{F[1],F[1]}-> {F[1],F[1]}, InsertionLevel -> {Classes},
		Model -> FileNameJoin[{"FermionEFT","FermionEFT"}],GenericModel->FileNameJoin[{"FermionEFT","FermionEFT"}]];
		
Paint[diags, ColumnsXRows -> {4, 1}, Numbering -> None,SheetHeader->None, ImageSize->{512,256}];

FCClearScalarProducts[]
SetMandelstam[s, t, u, p1, p2, -p3, -p4, 0, 0, 0, 0];
amp[0]=FCFAConvert[CreateFeynAmp[diags, Truncated -> True,PreFactor->1],
IncomingMomenta->{p1,p2},OutgoingMomenta->{p3,p4}, LoopMomenta->{q}, UndoChiralSplittings->True,
DropSumOver->True, List->False,SMP->True, ChangeDimension->D, Contract->True];

Could you please help resolve this issue?

[vsht]

You need to set Truncated to False, since the spinors are required to reconstruct the correct relative signs of the chains.

[vsht]

See also https://github.com/orgs/FeynCalc/discussions/221

[vsht]

FWIW I added a warning (1c30424) that should inform the user about this issue if the amplitude is truncated.

[Ichthyocentaur]

Thank you! Sorry for posting it in here, when it was not a bug but instead an error on my part.

[vsht]

No problem, the error message also was not particularly informative.

I'm closing this bug for now since all issues now appear to be resolved.