Double counting Fermion interchange with four Fermion interaction

Your Mathematica version

13.3.0 for Microsoft Windows (64-bit) (June 3, 2023)
Your FeynCalc version

Tested both 9.3.1 (hotfix branch) and 10.0.0 (master branch).
Did you try to reinstall FeynCalc (stable version) using the automatic installer to make sure that you have the latest bugfixes?

No, but both branches are up to date.
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.

(** User Mathematica initialization file **)
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.

The Feynman diagram from `FeynArts`:

FADiag = TopologyList[
  Process \[Rule] {F[
      3, {SequenceForm["Gen", 1], SequenceForm["Col", 1]}], -F[
       3, {SequenceForm["Gen", 2], SequenceForm[
        "Col", 2]}]} \[Rule] {F[2, {
SequenceForm["Gen", 3]}], -F[2, {
SequenceForm["Gen", 4]}]}, Model \[Rule] {"SMEFT-6-4F"}, 
  GenericModel \[Rule] {"SMEFT-6-4F"}, 
  InsertionLevel \[Rule] {Classes}, 
  ExcludeParticles \[Rule] {S, U, S[1], S[2], -S[3], S[3], -U[1], 
    U[1], -U[2], U[2], -U[4], U[4], -U[31], U[31], -U[32], 
    U[32], -U[4, {_}], U[4, {_}]}, ExcludeFieldPoints \[Rule] {}, 
  LastSelections \[Rule] {}][
 Topology[2][
   Propagator[Incoming][Vertex[1][1], Vertex[4][5], Field[1]], 
   Propagator[Incoming][Vertex[1][2], Vertex[4][5], Field[2]], 
   Propagator[Outgoing][Vertex[1][3], Vertex[3][6], Field[3]], 
   Propagator[Outgoing][Vertex[1][4], Vertex[3][6], Field[4]], 
   Propagator[Internal][Vertex[3][6], Vertex[3][7], Field[5]], 
   Propagator[FALoop[1]][Vertex[3][7], Vertex[4][5], Field[6]], 
   Propagator[FALoop[1]][Vertex[3][7], Vertex[4][5], 
    Field[7]]] \[Rule] 
  Insertions[Generic][
   FeynmanGraph[2, Generic \[Equal] 1][
     Field[1] \[Rule] 
      F[3, {SequenceForm["Gen", 1], SequenceForm["Col", 1]}], 
     Field[2] \[Rule] -F[
        3, {SequenceForm["Gen", 2], SequenceForm["Col", 2]}], 
     Field[3] \[Rule] -F[2, {
SequenceForm["Gen", 3]}], Field[4] \[Rule] F[2, {
SequenceForm["Gen", 4]}], Field[5] \[Rule] V, Field[6] \[Rule] F, 
     Field[7] \[Rule] F] \[Rule] 
    Insertions[Classes][
     FeynmanGraph[1, Classes \[Equal] 2][
      Field[1] \[Rule] 
       F[3, {SequenceForm["Gen", 1], SequenceForm["Col", 1]}], 
      Field[2] \[Rule] -F[
         3, {SequenceForm["Gen", 2], SequenceForm["Col", 2]}], 
      Field[3] \[Rule] -F[2, {
SequenceForm["Gen", 3]}], Field[4] \[Rule] F[2, {
SequenceForm["Gen", 4]}], Field[5] \[Rule] V[1], 
      Field[6] \[Rule] -F[
         3, {SequenceForm["Gen", 5], SequenceForm["Col", 5]}], 
      Field[7] \[Rule] 
       F[3, {SequenceForm["Gen", 5], SequenceForm["Col", 5]}]]]]]

The Feynman amplitude from `FeynArts`:

``` FAamp = FAFeynAmpList[ Process \[Rule] {{F[ 3, {SequenceForm["Gen", 1], SequenceForm["Col", 1]}], SequenceForm["p", 1], mudiag[SequenceForm["Gen", 1], SequenceForm["Col", 1]], {(2 Q)/ 3}}, {-F[3, {SequenceForm["Gen", 2], SequenceForm["Col", 2]}], SequenceForm["p", 2], mudiag[SequenceForm["Gen", 2], SequenceForm[ "Col", 2]], {-((2 Q)/3)}}} \[Rule] {{F[2, { SequenceForm["Gen", 3]}], SequenceForm["k", 1], mldiag[ SequenceForm["Gen", 3]], {-Q, LeptonNumber}}, {-F[2, { SequenceForm["Gen", 4]}], SequenceForm["k", 2], mldiag[ SequenceForm["Gen", 4]], {Q, -LeptonNumber}}}, Model \[Rule] {"SMEFT-6-4F"}, GenericModel \[Rule] {"SMEFT-6-4F"}, AmplitudeLevel \[Rule] {Classes}, ExcludeParticles \[Rule] {S, U, S[1], S[2], -S[3], S[3], -U[1], U[1], -U[2], U[2], -U[4], U[4], -U[31], U[31], -U[32], U[32], -U[4, {_}], U[4, {_}]}, ExcludeFieldPoints \[Rule] {}, LastSelections \[Rule] {}][ FAFeynAmp[ GraphID[Topology \[Equal] 1, Generic \[Equal] 1, Classes \[Equal] 1, Number \[Equal] 1], Integral[ SequenceForm["q", 1]], (1/(16 \[Pi]^4))FAFeynAmpDenominator[1/( SequenceForm["q", 1])^2, 1/(SequenceForm["q", 1] + SequenceForm["k", 1] + SequenceForm[ "k", 2])^2] "g"[SequenceForm["Lor", 1], SequenceForm[ "Lor", 2]] 1/(-( SequenceForm["k", 1]) - SequenceForm[ "k", 2])^2 (6 \[ImaginaryI] gc105 (Cuu[SequenceForm["Gen", 2], SequenceForm["Gen", 5], SequenceForm["Gen", 5], SequenceForm[ "Gen", 1]] + Cuu[SequenceForm["Gen", 5], SequenceForm["Gen", 1], SequenceForm["Gen", 2], SequenceForm["Gen", 5]]) IndexDelta[ SequenceForm["Col", 1], SequenceForm["Col", 5]] IndexDelta[ SequenceForm["Col", 2], SequenceForm["Col", 5]] IndexSum[ DiracObject["ga"[a], SequenceForm["om", Subscript["+"]]][SequenceForm["Dir", 2], SequenceForm[ "Dir", 8]] DiracObject["ga"[a], SequenceForm["om", Subscript["+"]]][SequenceForm["Dir", 6], SequenceForm["Dir", 1]], {a, 1, 4}] + 6 \[ImaginaryI] gc105 (Cuu[SequenceForm["Gen", 2], SequenceForm[ "Gen", 1], SequenceForm["Gen", 5], SequenceForm["Gen", 5]] + Cuu[SequenceForm["Gen", 5], SequenceForm["Gen", 5], SequenceForm["Gen", 2], SequenceForm["Gen", 1]]) IndexDelta[ SequenceForm["Col", 1], SequenceForm["Col", 2]] IndexSum[ DiracObject["ga"[a], SequenceForm["om", Subscript["+"]]][SequenceForm["Dir", 2], SequenceForm[ "Dir", 1]] DiracObject["ga"[a], SequenceForm["om", Subscript["+"]]][SequenceForm["Dir", 6], SequenceForm["Dir", 8]], {a, 1, 4}]) SumOver[SequenceForm["Col", 5], 3] SumOver[ SequenceForm["Gen", 5], 3] SumOver[SequenceForm["Col", 1], 3, External] SumOver[SequenceForm["Col", 2], 3, External] SumOver[ SequenceForm["Gen", 1], 3, External] SumOver[SequenceForm[ "Gen", 2], 3, External] SumOver[SequenceForm["Gen", 3], 3, External] SumOver[SequenceForm["Gen", 4], 3, External] DiracObject["gs"[-( SequenceForm["q", 1])]][SequenceForm["Dir", 5], SequenceForm[ "Dir", 6]] DiracObject[ "gs"[SequenceForm["q", 1] + SequenceForm["k", 1] + SequenceForm[ "k", 2]]][SequenceForm["Dir", 7], SequenceForm[ "Dir", 8]] DiracObject[FADiracSpinor[-( SequenceForm["p", 1]), 0]][ SequenceForm["Dir", 1]] DiracObject[FADiracSpinor[-( SequenceForm["p", 2]), 0]][ SequenceForm[ "Dir", 2]] DiracObject[FADiracSpinor[SequenceForm["k", 1], 0]][ SequenceForm[ "Dir", 3]] DiracObject[FADiracSpinor[SequenceForm["k", 2], 0]][ SequenceForm[ "Dir", 4]] (\[ImaginaryI] gc82 IndexDelta[SequenceForm["Gen", 3], SequenceForm["Gen", 4]] DiracObject["ga"[ SequenceForm["Lor", 1]], SequenceForm["om", Subscript["-"]]][SequenceForm["Dir", 3], SequenceForm[ "Dir", 4]] + \[ImaginaryI] gc82 IndexDelta[SequenceForm[ "Gen", 3], SequenceForm["Gen", 4]] DiracObject["ga"[ SequenceForm["Lor", 1]], SequenceForm["om", Subscript["+"]]][SequenceForm["Dir", 3], SequenceForm[ "Dir", 4]]) (\[ImaginaryI] gc83 DiracObject["ga"[ SequenceForm["Lor", 2]], SequenceForm["om", Subscript["-"]]][SequenceForm["Dir", 7], SequenceForm[ "Dir", 5]] + \[ImaginaryI] gc83 DiracObject["ga"[ SequenceForm["Lor", 2]], SequenceForm["om", Subscript["+"]]][SequenceForm["Dir", 7], SequenceForm["Dir", 5]])]] ```

FeynCalc amplitude:

FCamp = FCFAConvert[FAamp, IncomingMomenta -> {p1, p2}, 
  OutgoingMomenta -> {k1, k2}, LoopMomenta -> {l}, 
  ChangeDimension -> D, DropSumOver -> True, 
  UndoChiralSplittings -> True]

The output is

``` {(1/(8 \[Pi]^4)) 3 \[ImaginaryI] gc105 gc82 gc83 (Cuu[SequenceForm["Gen", 2], SequenceForm["Gen", 1], SequenceForm["Gen", 5], SequenceForm[ "Gen", 5]] + Cuu[SequenceForm["Gen", 5], SequenceForm["Gen", 5], SequenceForm[ "Gen", 2], SequenceForm["Gen", 1]]) DiracTrace[ GSD[k1 + k2 + l] . GA[6] . DiracGamma[LorentzIndex[Ind6088, D], D] . GSD[l] . GAD[Lor2]] Spinor[Momentum[k1, D], 0, 1] . GAD[Lor1] . Spinor[-Momentum[k2, D], 0, 1] Spinor[-Momentum[p2, D], 0, 1] . DiracGamma[LorentzIndex[Ind6088, D], D] . GA[6] . Spinor[Momentum[p1, D], 0, 1] FAD[-k1 - k2] FAD[l, k1 + k2 + l] IndexDelta[SequenceForm["Gen", 3], SequenceForm[ "Gen", 4]] MTD[Lor1, Lor2] SDF[Col1, Col2] + (1/(8 \[Pi]^4)) 3 \[ImaginaryI] gc105 gc82 gc83 (Cuu[SequenceForm["Gen", 2], SequenceForm["Gen", 5], SequenceForm["Gen", 5], SequenceForm[ "Gen", 1]] + Cuu[SequenceForm["Gen", 5], SequenceForm["Gen", 1], SequenceForm[ "Gen", 2], SequenceForm["Gen", 5]]) Spinor[Momentum[k1, D], 0, 1] . GAD[Lor1] . Spinor[-Momentum[k2, D], 0, 1] Spinor[-Momentum[p2, D], 0, 1] . (-DiracGamma[LorentzIndex[Ind6087, D], D] . GA[6] . GSD[k1 + k2 + l]) . GAD[Lor2] . (-GSD[l]) . DiracGamma[LorentzIndex[Ind6087, D], D] . GA[6] . Spinor[Momentum[p1, D], 0, 1] FAD[-k1 - k2] FAD[l, k1 + k2 + l] IndexDelta[SequenceForm["Gen", 3], SequenceForm[ "Gen", 4]] MTD[Lor1, Lor2] SDF[Col1, Col5] SDF[Col2, Col5]} ```

For this diagram, the Feynman rule is: If one naively inserts this Feynman rule into the diagram, the result is exactly what FeynCalc produces. But when deriving this Feynman rule, the interchange of 4 identical Fermion fields is considered, hence the minus sign in the second term. This interchange is also accounted for when producing the fermion loop in the diagram above. Therefore, one should not double-count the minus sign when a four-fermion interaction with four identical fermions is involved.

The correct form of the amplitude is similar to FCamp, but the term without trace should have an additional minus sign.

In my case applying FCFAConvert to FAamp yields

{("g"["Lor"1, "Lor"2]*FeynAmpDenominator["q"1^(-2), ("k"1 + "k"2 + "q"1)^(-2)]*DiracObject["gs"[-"q"1]]["Dir"5, "Dir"6]*DiracObject["gs"["k"1 + "k"2 + "q"1]]["Dir"7, "Dir"8]*DiracObject[Spinor[Momentum["k"1, D], 0, 1]]["Dir"3]*
   DiracObject[Spinor[Momentum["k"2, D], 0, 1]]["Dir"4]*DiracObject[Spinor[-Momentum["p"1, D], 0, 1]]["Dir"1]*DiracObject[Spinor[-Momentum["p"2, D], 0, 1]]["Dir"2]*
   ((6*I)*gc105*(Cuu["Gen"2, "Gen"5, "Gen"5, "Gen"1] + Cuu["Gen"5, "Gen"1, "Gen"2, "Gen"5])*IndexDelta["Col"1, "Col"5]*IndexDelta["Col"2, "Col"5]*DiracObject["ga"[Ind19], "om"Subscript["+"]]["Dir"2, "Dir"8]*
     DiracObject["ga"[Ind19], "om"Subscript["+"]]["Dir"6, "Dir"1] + (6*I)*gc105*(Cuu["Gen"2, "Gen"1, "Gen"5, "Gen"5] + Cuu["Gen"5, "Gen"5, "Gen"2, "Gen"1])*IndexDelta["Col"1, "Col"2]*DiracObject["ga"[Ind20], "om"Subscript["+"]]["Dir"2, "Dir"1]*
     DiracObject["ga"[Ind20], "om"Subscript["+"]]["Dir"6, "Dir"8])*(I*gc82*IndexDelta["Gen"3, "Gen"4]*DiracObject["ga"["Lor"1], "om"Subscript["-"]]["Dir"3, "Dir"4] + 
    I*gc82*IndexDelta["Gen"3, "Gen"4]*DiracObject["ga"["Lor"1], "om"Subscript["+"]]["Dir"3, "Dir"4])*(I*gc83*DiracObject["ga"["Lor"2], "om"Subscript["-"]]["Dir"7, "Dir"5] + I*gc83*DiracObject["ga"["Lor"2], "om"Subscript["+"]]["Dir"7, "Dir"5]))/
  (16*Pi^4*(-"k"1 - "k"2)^2)}

so there must be something wrong with the input.

Apart from that, you can try to create a model where two fermion lines are connected via a fake particle and check what comes out. This is how 4-fermion vertices are normally treated.

Then you can also cross check this with FormCalc, which implements the same 4-fermion resolving algorithm that FeynCalc uses.

For the moment I don't really see an issue here.

I think there's an issue with copying directly the StandardForm output with CopyToClipboard. Compress it should work:

FAamp = Uncompress@"1:eJzlWFtv0zAUbpN2K2NM04TgceMHTOIuXkcv06SWlbX9AV7qZoY0LrmMjt+OAF+W2HHdxOnCNEQfrMbnfMfn7pO8uMQXs3qtVpttkyV8QpbeSQ/e+CfzRR+F0cyi2w2yXMQeDCnTMMAODMMsRTwxVCPzxGj0kB6y6VHWKrVJljN/CpfhDvl3Cn0YgAhhH9VTgGDZIv/a2MNxIJF3qXyyNcBz6EfxPGwxfgfPke9KfBQ8j6cIuBWcXNdbMkZzeOsgqsUFEwg8RHeYC5g3PuvclMLRb/K7k+usfAMEOd91gs/YdSYnl3TdL+IN2XWG+ZZ3lFZ14WWtV87jyMVpQtVTr3gGXskLoAg5t4+d3oeLCPuf4vklDMrlSmmzG4ZmW+XNbqyYLYdQMX7F7Gz/odwDPIVe1rgRPXI06PbGx++P3/YU0G6iD3JKY/fIH9IKPRTFU9iH1yqadcS2B8IQKh0x3Cd/ukvHI8AhCCLkkL3ZwYofRmydcLFsQ8ot/pwtJyniEo+tYOxCzEQ5Z6J0PD3GUjDFuk14BqjPBZhDBXNogDlSMEfrMJZOt7RSagniowf8rxswpylwIFKgh6A3HZIaikKOyRZmJuH6gKQl9KBD60dlZz3wUU26p3lroIl4GoDF1VlH2Nz9FgOPVfAYL0jzdW+UIHOGbVEha+m3Sa6l017A61XKKX6JRNANgJd7z0QwoB2eIZ9loqW5BygX2s/oMMTfYRDSxyFCP9P4PpV91IE+vc1AhG9bynNGJbPMArh0d4WhWFlUS9XIk2Yn4Rt6cWjsCNMryLBnS7ruMV0HMCLRHkM/TAwWDZyGu48DIuqHdgqS6ZaRD6xVH+hK0tTe8ljFB6kqW5p8Y9mOSdOHSwZBKZPrvHr5TpHAio/GuB3HG4xFCkuzEhbmqRKKmYqs0LydhKUDvQgYTfha8gYSC6ZTIbGVkEfxXDeh8r3HZOmgADjnl19Iy+bNj79Ise0BIJW2ZDcI4ETemdpXhEoGixtSMhSI9bY2E/FavQW1lVDvWautXCrDph0wHXoAb6DsPm0+sBI0yHTTKq22BA1YSplXVQmKQ/+pgkmxD6pgWsUFk1YI8fL5NVy5vpVeJr0K5IA0ucaBtvlp9QTEEqG75PNMGQnWZhI0KX53KdXoYlcipaGXIl7GtNm7K7J3RKb3q5LDk2Yg1WWsri1psr2Mng9iZF5v03ZxBdfXWsxwUl8ZLZAvj+CGoZG+a0pz7do+d1/6WMX6WJvrUxzcgrPtv3W2gd3Z72BicCkafdLvDa7z4bXJzb32s2PxJ7oyV2HhW2P+JSkyam2g8l35PzsvvyMbOC+bh3ZZV76p1GCr2OD8bMlvyOyCevgm5se02MQ/vAEcYw=="

Can you try again with this definition?

I don't think a fake intermediate particle can correctly produce the fermion chain in this case. I'll look into FormCalc as suggested.

Thanks, now it works. With

amp = FCFAConvert[FAamp, IncomingMomenta -> {p1, p2}, 
   OutgoingMomenta -> {k1, k2}, LoopMomenta -> {l}, 
   ChangeDimension -> D, DropSumOver -> True, 
   UndoChiralSplittings -> True, FCFADiracChainJoin -> False] // FCE

you get the amplitude with explicit Dirac indices, where the relative signs are still unfixed and some indices are reversed.

{-1/16*(gc82*gc83*DCHN[Spinor[Momentum[k1, D], 0, 1], Dir3]*DCHN[Spinor[Momentum[k2, D], 0, 1], Dir4]*DCHN[Spinor[-Momentum[p1, D], 0, 1], Dir1]*DCHN[Spinor[-Momentum[p2, D], 0, 1], Dir2]*
    DCHN[GAD[Lor1], Dir3, Dir4]*DCHN[GAD[Lor2], Dir7, Dir5]*DCHN[-GSD[l], Dir5, Dir6]*DCHN[GSD[k1 + k2 + l], Dir7, Dir8]*FAD[-k1 - k2]*FAD[l, k1 + k2 + l]*IndexDelta[Index[Generation, 3], Index[Generation, 4]]*
    MTD[Lor1, Lor2]*((6*I)*gc105*(Cuu[Index[Generation, 2], Index[Generation, 1], Index[Generation, 5], Index[Generation, 5]] + Cuu[Index[Generation, 5], Index[Generation, 5], Index[Generation, 2], 
        Index[Generation, 1]])*DiracChain[DiracGamma[LorentzIndex[Ind194, D], D] . GA[6], DiracIndex[Dir2], DiracIndex[Dir1]]*DiracChain[DiracGamma[LorentzIndex[Ind194, D], D] . GA[6], DiracIndex[Dir6], 
       DiracIndex[Dir8]]*SDF[Col1, Col2] + (6*I)*gc105*(Cuu[Index[Generation, 2], Index[Generation, 5], Index[Generation, 5], Index[Generation, 1]] + 
       Cuu[Index[Generation, 5], Index[Generation, 1], Index[Generation, 2], Index[Generation, 5]])*DiracChain[DiracGamma[LorentzIndex[Ind193, D], D] . GA[6], DiracIndex[Dir2], DiracIndex[Dir8]]*
      DiracChain[DiracGamma[LorentzIndex[Ind193, D], D] . GA[6], DiracIndex[Dir6], DiracIndex[Dir1]]*SDF[Col1, Col5]*SDF[Col2, Col5]))/Pi^4}

Then FCFADiracChainJoin overtakes by building up closed spinor chains and fixing relative signs

https://github.com/FeynCalc/feyncalc/blob/master/FeynCalc/Feynman/FCFADiracChainJoin.m

FCFADiracChainJoin[amp, FCE -> True]
(*
{(((3*I)/8)*gc105*gc82*gc83*(Cuu[Index[Generation, 2], Index[Generation, 1], Index[Generation, 5], Index[Generation, 5]] + Cuu[Index[Generation, 5], Index[Generation, 5], Index[Generation, 2], 
      Index[Generation, 1]])*DiracTrace[GSD[k1 + k2 + l] . GA[6] . GAD[Ind232] . GSD[l] . GAD[Lor2]]*Spinor[Momentum[k1, D], 0, 1] . GAD[Lor1] . Spinor[-Momentum[k2, D], 0, 1]*
    Spinor[-Momentum[p2, D], 0, 1] . GAD[Ind232] . GA[6] . Spinor[Momentum[p1, D], 0, 1]*FAD[-k1 - k2]*FAD[l, k1 + k2 + l]*IndexDelta[Index[Generation, 3], Index[Generation, 4]]*MTD[Lor1, Lor2]*SDF[Col1, Col2])/
   Pi^4 + (((3*I)/8)*gc105*gc82*gc83*(Cuu[Index[Generation, 2], Index[Generation, 5], Index[Generation, 5], Index[Generation, 1]] + Cuu[Index[Generation, 5], Index[Generation, 1], Index[Generation, 2], 
      Index[Generation, 5]])*Spinor[Momentum[k1, D], 0, 1] . GAD[Lor1] . Spinor[-Momentum[k2, D], 0, 1]*Spinor[-Momentum[p2, D], 0, 1] . (-GAD[Ind231] . GA[6] . GSD[k1 + k2 + l]) . GAD[Lor2] . (-GSD[l]) . 
     GAD[Ind231] . GA[6] . Spinor[Momentum[p1, D], 0, 1]*FAD[-k1 - k2]*FAD[l, k1 + k2 + l]*IndexDelta[Index[Generation, 3], Index[Generation, 4]]*MTD[Lor1, Lor2]*SDF[Col1, Col5]*SDF[Col2, Col5])/Pi^4}*)

It uses the same algorithm as FormCalc, so the final results should be identical. And FormCalc is really trustworthy in that respect. Unless your FeynArts model itself has an issue in it.

I don't think a fake intermediate particle can correctly produce the fermion chain in this case.

I'm pretty sure it can because this is how people handle this complication for decades in QGRAF, MadGraph etc. A fake particle is actually the most clever approach because it gives you correct relative signs for free.

Apology for the late reply. The result you are getting is exactly where I claim the relative signs are wrong. There are two terms in the amplitude: one with trace and one without. There should be an extra minus relative sign between them. However, I'm not exactly sure what's the problem.

Also, thanks for the suggestion on FormCalc. FormCalc actually produces correct signs while FeynCalc cannot. And you are absolutely correct about the heavy mediator. Initially, I forgot the t-channel diagrams that give the 2nd term in my Feynman rules after integrating out heavy modes, which is why I had doubts.

One other thing I noticed is the missing minus sign from a closed fermion loop. Did you implement those in FCFADiracChainJoin? I tested a little bit with the SM model shipped with FeynArts and my model that uses DiracChain. Regular diagrams such as a vertex correction produce the same sign, but a closed fermion loop (i.e. a bubble) gives different signs for two models.

If you can provide the full FA model and some example I could look into it in more details. FeynCalc should actually give the same results as FormCalc, since the underlying algorithm is the same.

FeynCalc / feyncalc

Double counting Fermion interchange with four Fermion interaction #234