Issue with γ-matrix simplifications after FermionSpinSum in FeynCalc

[impermast]

Hello,

I am experiencing an issue with FeynCalc (version 9.3) in Mathematica while computing the amplitude for the process tt → Z → ZZ in the 2HDM model.

Previously, I successfully computed the amplitude for ZZZ, and γ-matrix expressions were correctly simplified. However, after introducing quarks in the initial state and applying FermionSpinSum, along with SUNSimplify, γ-matrix expressions no longer simplify as expected.

MWE

<< FeynArts`
<< FeynCalc`

(*only one diagram for simplicity*)
diagsAll = InsertFields[
    CreateTopologies[1, 2 -> 2, ExcludeTopologies -> {Tadpoles}], 
    {F[9], -F[9]} -> {V[2], V[2]}, 
    InsertionLevel -> {Classes}, 
    Model -> THDMCPV, 
    ExcludeParticles -> {F[_], V[1], V[3], S[2], S[3], U[_]}, 
    LastSelections -> {V[2]}
];

diagHHH = DiagramExtract[diagsAll, {19}];

ampHHH[0] = FCFAConvert[
    CreateFeynAmp[diagHHH, Truncated -> True], 
    IncomingMomenta -> {k1, k2}, 
    OutgoingMomenta -> {p1, p2}, 
    LoopMomenta -> {q}, 
    TransversePolarizationVectors -> {p1, p2}, 
    UndoChiralSplittings -> True, 
    ChangeDimension -> D, 
    List -> False, 
    SMP -> True, 
    Contract -> False, 
    DropSumOver -> True
];

ampHHH[1] = ampHHH[0] // FCTraceFactor // SUNSimplify;

FCClearScalarProducts[];
SetMandelstam[s, t, u, k1, k2, -p1, -p2, SMP["m_t"], SMP["m_t"], SMP["m_Z"], SMP["m_Z"]];

ampHHH[2] = ampHHH[1] // TID[#, q, ToPaVe -> True, UsePaVeBasis -> True] &;

ampHHH[3] = (ampHHH[2] /. M$FACouplings) // FCReplaceD[#, D -> 4 - 2 Epsilon] & //
    Series[#, {Epsilon, 0, 0}] & // Normal;

ampHHH[4] = ampHHH[3] // ChangeDimension[#, 4] &;

(*just contract amp, no summing*)
(ampHHH[4] ComplexConjugate[ampHHH[4]]) // Contract

(*fermion and vector sum*)
ampSquared[0] = 1/2 (ampHHH[4] (ComplexConjugate[ampHHH[4]])) //  
                FermionSpinSum // Contract//
        DiracSimplify //
        SUNSimplify  // 
        DoPolarizationSums[#, k1] & // 
        DoPolarizationSums[#, k2] &//
        TrickMandelstam[#, {s, t, u, 2 SMP["m_t"]^2 + 2 SMP["m_Z"]^2}] &;

ampSquared[0]

Expected Behavior

The squared amplitude should not explicitly depend on γ-matrices. After applying the necessary simplifications, terms involving DiracGamma should be completely eliminated, leaving a scalar result.

Observed Behavior

In ampSquared[0], I still have terms of the form:

DiracGamma[Momentum[p1]] \cdot DiracGamma[7] \cdot DiracGamma[Momentum[p2]]

These terms exist even after applying FermionSpinSum, DiracSimplify, SUNSimplify, and other simplifications.

Questions:

1. What is the correct order of operations to ensure proper γ-matrix simplification?

• Contract
• TID (Tensor Integral Decomposition)
• FermionSpinSum
• DiracSimplify
• SUNSimplify
• DoPolarizationSums[#, k1] & // DoPolarizationSums[#, k2] & //

2. Should DiracSimplify be applied before or after FermionSpinSum to ensure correct simplifications?

I would be grateful for any clarification on whether I am missing a step or misapplying a function. Thank you in advance for your help!