FeynCalc / feyncalc

Mathematica package for algebraic calculations in elementary particle physics.
https://feyncalc.github.io
GNU General Public License v3.0
310 stars 87 forks source link

Bug in FCMultiLoopTID: Mismatch between pure-evanescent or 4D structures, vs. fully d-D structures #52

Closed HBelusca closed 4 years ago

HBelusca commented 4 years ago

11.3.0 for Microsoft Windows (64-bit) (March 7, 2018)

9.4.0 master (a3aab44faccb7e68d0fab29a626b1ddc47667725)

Using git master repo checkout


Using FCMultiLoopTID[] returns unexpected results when the provided amplitude contains explicit evanescent (in D-4 dimensions) or 4-dimensional structures. This causes problems in analytical calculation of (2+) loop amplitudes in e.g. BMHV scheme. This problem also shows up, at 1-loop level, when comparing the tensor reduction obtained via the usual TID[] function, and by FCMultiLoopTID[] with 1 momentum specification.

The observed effect is that it appears, that somehow FCMultiLoopTID[] "replaces" all occurrences of loop momenta (when evanescent or 4-dimensional projections of these are present) to their fully d-dimensional counterparts. See the attached notebook (in zip archive) for explicit examples, as well as for a "workaround".

FeynCalc_bug_FCMultiLoopTID.zip

HBelusca commented 4 years ago

(Note: The "FCMultiLoopTID::nomulti" warning in LoopIntegrals/FCMultiLoopTID.m does not seem to be used anywhere in FeynCalc.)

HBelusca commented 4 years ago

After much observation in the code (file: LoopIntegrals/FCMultiLoopTID.m), the bug seems to be in the code block starting with

FCPrint[1, "FCMultiLoopTID: Handling 4 and D-4 dimensional loop momenta.", FCDoControl->mltidVerbose];

and/or what the amplitude form expectation when the following uncontraction step:

ex = Uncontract[ex, Sequence@@qs, Pair -> optUncontract, CartesianPair-> optUncontract, FCI->True];

is done few lines above. Indeed I traced the problem to the fact that scalar products involving loop momenta in either 4 or 4-D dimensions can as well be present in the amplitude, and those are not handled at all either by this Uncontract step, and/or by the step "Handling 4 and D-4 dimensional loop momenta" just after.

As an alternative to the workaround done in the notebook attached in the zip file, I've found another simpler one:

FCMultiLoopTID[amplitude, {loop, momenta}, Uncontract -> {loop, momenta}]

that again, might be limited to the fact that now a lot of uncontracted structures might appear, and cause trouble for the helper functions Tdec (e.g. obtaining the error "Tdec::indices: Tdec has detected that the integral contains loop momenta with idential (sic.) Lorentz indices. Evaluation aborted!"), as well as a dramatic performance drop.

vsht commented 4 years ago

Thanks for the very interesting bug report. I had to think a bit why this happens, but after all it is clear why the function fails.

FCMultiLoopTID is never trying to perform a full reduction in the style of TID. The reason is that this makes no sense beyond 1-loop. Multiloop integrals may have irreducible numerators that can be removed (or rather traded for higher denominator powers) only via IBP reduction. Tensor reduction won't help here.

So FCMultiLoopTID only performs reduction w.r.t momenta with free indices or those contracted to polarization vectors, eps tensors or Dirac matrices (or those you specify with the Uncontract option, but this is rather reserved for special cases). That is, all kind of stuff that doesn't belong into a scalar loop integral. The output is then in the proper form for the subsequent steps in a multiloop calculation: topology identification and IBP reduction.

As far as the reduction of integrals with mixed dimensions is concerned, what FeynCalc does here is the following:

4-dim loop momenta are rewritten as \bar{l}^{mu} = \bar{g}^{mu nu} l_{nu} (D-4)-dim loop momenta are rewritten as \hat{l}^{mu} = \hat{g}^{mu nu} l_{nu}

So in both cases the loop momentum is temporarily "upgraded" to D dimensions. The external momenta in the propagators are assumed to be D-dimensional as well. This way you can recast a mixed integral into a pure D-dimensional integral and then tensor reduce it the usual way.

Now the issue with FCMultiLoopTID is that it doesn't apply this procedure to all loop momenta in the numerator but only to those that it selects according to the criteria I described above. However, in this case the integral you are reducing is not purely D-dimensional! We still may have SPE- or SP-type scalar products (as in your example) and thus the whole procedure becomes ill-defined. Notice that this doesn't happen to TID, as it always uncontract all loop momenta in the numerator.

The fix is essentially what you have already observed in your code: In the case of mixed integrals FCMultiLoopTID should also uncontract all loop momenta in the numerators, even those contracted to other momenta. Otherwise the output will not be correct.

I've pushed a fix that should handle this issue while still preserving the original behavior of FCMultiLoopTID not to touch D-dim scalar products. Let me know if it works for you.

As for the last section of your notebook, well sure SPD[q1, q2] FVE[q2, mu] FAD[q1, q1 - p, q2] is zero, it is after all anitsymmetric in q2, right?

HBelusca commented 4 years ago

Thank you very much! Your fix works perfectly, therefore the report can be closed.

My alternative I was experimenting with, was to do the following:

diff --git a/FeynCalc/LoopIntegrals/FCMultiLoopTID.m b/FeynCalc/LoopIntegrals/FCMultiLoopTID.m
index a4df906c..d90d6259 100644
--- a/FeynCalc/LoopIntegrals/FCMultiLoopTID.m
+++ b/FeynCalc/LoopIntegrals/FCMultiLoopTID.m
@@ -173,22 +174,25 @@ FCMultiLoopTID[expr_/;Head[expr]=!=List, qs_List/; FreeQ[qs, OptionQ], OptionsPa
                FCPrint[1, "FCMultiLoopTID: Done uncontracting Lorentz indices, timing: ", N[AbsoluteTime[] - time, 4], FCDoControl->mltidVerbose];
                FCPrint[3, "FCMultiLoopTID: After Uncontract: ", ex, FCDoControl->mltidVerbose];

-               If[ (FeynCalc`Package`DiracGammaScheme === "BMHV") && !FreeQ2[ex,{LorentzIndex,CartesianIndex}],
+               If[ (FeynCalc`Package`DiracGammaScheme === "BMHV") && !FreeQ2[ex,{LorentzIndex,CartesianIndex,Momentum}],
                        time=AbsoluteTime[];
                        FCPrint[1, "FCMultiLoopTID: Handling 4 and D-4 dimensional loop momenta.", FCDoControl->mltidVerbose];
-                       ex = ex /. {
-                               Pair[Momentum[q_,n-4],LorentzIndex[i_,n-4]]/; MemberQ[qs,q]:>
-                                       (tmpli=Unique[];  Pair[Momentum[q,n],LorentzIndex[tmpli,n]] Pair[LorentzIndex[tmpli,n-4],LorentzIndex[i,n-4]]),
-                               Pair[Momentum[q_],LorentzIndex[i_]]/; MemberQ[qs,q] :>
-                                       (tmpli=Unique[];  Pair[Momentum[q,n],LorentzIndex[tmpli,n]] Pair[LorentzIndex[tmpli],LorentzIndex[i]]),
-
-                               CartesianPair[CartesianMomentum[q_,n-4],CartesianIndex[i_,n-4]]/; MemberQ[qs,q] :>
-                                       (tmpli=Unique[];  CartesianPair[CartesianMomentum[q,n-1],CartesianIndex[tmpli,n-1]] CartesianPair[CartesianIndex[tmpli,n-4],CartesianIndex[i,n-4]]),
-                               CartesianPair[CartesianMomentum[q_],CartesianIndex[i_]]/; MemberQ[qs,q] :>
-                                       (tmpli=Unique[];  CartesianPair[CartesianMomentum[q,n-1],CartesianIndex[tmpli,n-1]] CartesianPair[CartesianIndex[tmpli],CartesianIndex[i]])
+                       ex = ex //. {
+                               Pair[Momentum[q_,n-4],r_] /; MemberQ[qs,q] :>
+                                       (tmpli=Unique[];  Pair[Momentum[q,n],LorentzIndex[tmpli,n]] Pair[LorentzIndex[tmpli,n-4],r]),
+                               Pair[Momentum[q_],r_] /; MemberQ[qs,q] :>
+                                       (tmpli=Unique[];  Pair[Momentum[q,n],LorentzIndex[tmpli,n]] Pair[LorentzIndex[tmpli],r]),
+
+                               CartesianPair[CartesianMomentum[q_,n-4],r_] /; MemberQ[qs,q] :>
+                                       (tmpli=Unique[];  CartesianPair[CartesianMomentum[q,n-1],CartesianIndex[tmpli,n-1]] CartesianPair[CartesianIndex[tmpli,n-4],r]),
+                               CartesianPair[CartesianMomentum[q_],r_] /; MemberQ[qs,q] :>
+                                       (tmpli=Unique[];  CartesianPair[CartesianMomentum[q,n-1],CartesianIndex[tmpli,n-1]] CartesianPair[CartesianIndex[tmpli],r])
                        };

which ultimately does the same thing (for e.g. evanescent scalar product between two loop momenta), but would not introduce extra (e.g. evanescent) metric tensors in the case of a (e.g. evanescent) scalar product between a loop momentum and an external momentum: Pair[Momentum[k1,D-4], Momentum[q1,D-4]] would give Pair[Momentum[k1,D-4], LorentzIndex[mu,D-4]] Pair[Momentum[q1,D-4], LorentzIndex[mu,D-4]] only, instead of the longer (but equivalent) Pair[Momentum[k1,D-4], LorentzIndex[mu,D-4]] Pair[Momentum[q1,D], LorentzIndex[rho,D]] Pair[LorentzIndex[rho,D-4], LorentzIndex[mu,D-4]] as what now happens.

As for the last section of your notebook, well sure SPD[q1, q2] FVE[q2, mu] FAD[q1, q1 - p, q2] is zero, it is after all anitsymmetric in q2, right?

You are correct!

HBelusca commented 4 years ago

By the way (unrelated to this bug at hand), in the FCMultiLoopTID function, are these momenta called "q" correct, in (~ line 130):

(h: Pair|FeynAmpDenominator)[x__] /; !FreeQ[{x}, q] :> FRH[h[x], IsolateNames->mltidIsolate];

and in (~ lines 188-189)

            If[ !FreeQ2[ex, {Pair[Momentum[q,n-4],LorentzIndex[_,n-4]],Pair[Momentum[q],LorentzIndex[_]],
                    CartesianPair[CartesianMomentum[q,n-4],CartesianIndex[_,n-4]],CartesianPair[CartesianMomentum[q],CartesianIndex[_]]}],

? Shouldn't they be instead the loop momenta (e.g. Alternatives@@qs for the latter case) ?

vsht commented 4 years ago

Should be also taken care of. The alternative you were experimenting with can produce tensor integrals with identical indices, as in the case of SPE[q1,q1]. Actually, now I see why I had the corresponding warning in Tdec

Things like Tdec[{{l, mu1}, {l, mu1}, {l, mu3}}, {p}] do not work correctly. It is much simpler to calculate Tdec[{{l, mu1}, {l, mu2}, {l, mu3}}, {p}] and contract withg^{mu nu} afterwards.

vsht commented 4 years ago

Since there seem to be no further issues related to the fix, I'm closing this bug.