Wrong result for three-dimensional loop integral with a single term in the F-polynomial

[apik]

In the attached example pysecdec-se2l3d.zip trying to calculate two-loop massless sunset in d=3-2*ep, which is logarithmically divergent, its pole part is dependent on external momenta squared. If I rewrite the same example using auxiliary mass at the preparation step and set m=0 at the integration step the result is correct J=-Pi/eps+O(eps^0).

m=0 from the beginning:

J(-0.1) =  + (-3.14159359433898888e-01 +/- 2.66294995075465029e-05)*eps^-1 + O(eps^0)
J(-0.5) =  + (-1.57079679716949472e+00 +/- 1.33147497538035772e-04)*eps^-1 + O(eps^0)
J(-1.0) =  + (-3.14159359433898944e+00 +/- 2.66294995076071545e-04)*eps^-1 + O(eps^0)
J(-1.5) =  + (-4.71239039150848349e+00 +/- 3.99442492613653010e-04)*eps^-1 + O(eps^0)
J(-2.0) =  + (-6.28318718867797887e+00 +/- 5.32589990152143090e-04)*eps^-1 + O(eps^0)

m=0 during integration:

J(-0.1) =  + (-3.14159359433898855e+00 +/- 2.66294995075954343e-04)*eps^-1 + O(eps^0)
J(-0.5) =  + (-3.14159359433898855e+00 +/- 2.66294995075954343e-04)*eps^-1 + O(eps^0)
J(-1.0) =  + (-3.14159359433898855e+00 +/- 2.66294995075954343e-04)*eps^-1 + O(eps^0)
J(-1.5) =  + (-3.14159359433898855e+00 +/- 2.66294995075954343e-04)*eps^-1 + O(eps^0)
J(-2.0) =  + (-3.14159359433898855e+00 +/- 2.66294995075954343e-04)*eps^-1 + O(eps^0)

The only difference between these two considered examples is the form of the F-polynomial, which has a single term in the first case.

[magv]

Hi, Andrey. It seems that if you'll enable contour_deformation = True, then you get -Pi/eps in all cases. Not sure what's going on without it; we'll look into it.

[apik]

I'm interested in the calculation of the Euclidean kinematics integral with F>=0, so contour_deformation = True not seems to be an appropriate solution due to performance reasons.

The same story with the 3-loop sunset in its logarithmic dimension d=8/3-2*ep see attached example pysecdec-se3l83d.zip, correct result is Gamma[1/3]^3/eps+O(eps^0)

m=0 from the beginning:

J(-0.1) =  + (1.92457620986330835e+00 +/- 1.21831321541757544e-03)*eps^-1 + O(eps^0)
J(-0.5) =  + (9.62288104931645094e+00 +/- 6.09156607702509602e-03)*eps^-1 + O(eps^0)
J(-1.0) =  + (1.92457620986329019e+01 +/- 1.21831321540501920e-02)*eps^-1 + O(eps^0)
J(-1.5) =  + (2.88686431479495447e+01 +/- 1.82746982312152577e-02)*eps^-1 + O(eps^0)
J(-2.0) =  + (3.84915241972658038e+01 +/- 2.43662643081003841e-02)*eps^-1 + O(eps^0)

m!=0:

J(-0.1) =  + (1.92457620986329019e+01 +/- 1.21831321540501920e-02)*eps^-1 + O(eps^0)
J(-0.5) =  + (1.92457620986329019e+01 +/- 1.21831321540501920e-02)*eps^-1 + O(eps^0)
J(-1.0) =  + (1.92457620986329019e+01 +/- 1.21831321540501920e-02)*eps^-1 + O(eps^0)
J(-1.5) =  + (1.92457620986329019e+01 +/- 1.21831321540501920e-02)*eps^-1 + O(eps^0)
J(-2.0) =  + (1.92457620986329019e+01 +/- 1.21831321540501920e-02)*eps^-1 + O(eps^0)

[magv]

No, I'm not yet sure if contour_deformation = True is the solution. Just investigating.

---

In the first case pySecDec tries to integrate

(+(6))*((U)^(+(-3/2) + (3)*eps))*(+(-s))/Sqrt[x2]

while in the second case it is

(+(6))*((U)^(+(-3/2) + (3)*eps))*((F)^(+(-2)*eps))/Sqrt[x1]

where U is 1 + x2 + x1*x2 and F is -s*x1*x2 (ignoring terms proportional to m).

This is how you get the s dependence in the first case and -Pi/eps in the second.

[spj101]

Thanks @apik for reporting this. I have pushed a fix to the proposed branch. I'd be happy if you confirm that this works for you (also on any other examples you may have).

Apparently, there was a long-standing bug that caused exponents in the original integrand expression to be modified when derivatives were computed. This bug is only triggered by certain simple expressions but I don't yet understand under which conditions it occurs.

In any case, clearly, the original integrand should not be modified and it is sensible (and safer) to take a copy when computing derivatives. This could perhaps have a performance and memory impact for huge expressions. We may profile this and adjust our fix in future versions.