search

gudrunhe / secdec

a program to evaluate dimensionally regulated parameter integrals numerically
GNU General Public License v3.0
25 stars 5 forks source link

Issues compiling a 3-line 2-loop eikonal integral #27

Open vsht opened 2 years ago

vsht commented 2 years ago

Hi,

I'd like to use pySecDec to check my analytic result for the following on-shell eikonal integral with $q=m_1^2$ and $m_1>m_2$

Bildschirmfoto vom 2022-08-01 12-33-20

The $i \eta$-prescription automatically gets rewritten to match the pySecDec-convention.

My generation script looks as follows:

#!/usr/bin/env python3
from pySecDec import sum_package, loop_package, loop_regions, LoopIntegralFromPropagators
import pySecDec as psd

li = LoopIntegralFromPropagators(
['-m2**2 + p1**2', '2*p3*q', '-m2**2 + (-p1 + p3)**2'],
loop_momenta = ['p1', 'p3'],
powerlist = [1, 1, 1],
dimensionality = '4 - 2*eps',
Feynman_parameters = 'x',
replacement_rules = [('q**2','m1**2')],
regulators = ['eps']
)

loop_package(
name = 'loopint',
loop_integral = li,
requested_orders = [2],
real_parameters = ['m1', 'm2'],
contour_deformation = True,
additional_prefactor = '(-1)*exp(2*EulerGamma*eps)',
decomposition_method = 'geometric',
form_work_space = '200M'
)

If I omit form_work_space, then the compilation stops progressing quite soon and upon hitting Ctrl-C I see 

=== Workspace overflow. 50000000 bytes is not enough.
=== Change parameter WorkSpace in form.set

However, even with a larger workspace the compilation still doesn't work and after 15-20 minutes I end up with 

=== Workspace overflow. 200000000 bytes is not enough.

I think it is related to the contour deformation, since I've already observed a similar behavior with other eikonal integrals. However, in this computation I know that some of my masters will have an imaginary part and so I cannot disable the contour deformation for all integrals.

Is there some trick to get the compilation through?

FYI, FIESTA 5 calculates this integral within seconds:

SetOptions[FIESTA, "NumberOfSubkernels" -> 4, "NumberOfLinks" -> 4, 
  "ComplexMode" -> True, "ReturnErrorWithBrackets" -> True, 
  "Integrator" -> "quasiMonteCarlo", "Precision" -> 14(*,
  "ZeroCheckCount" -> 1000*)];

SDEvaluate[
 UF[{p1, p3}, {m2^2 - p1^2, -2*p3*q, -(p3 - p1)^2 + m2^2}, {q^2 -> 
    m1^2, m1 -> 2.76, m2 -> 1.421}], {1, 1, 1}, 2]

(*ep^2 (1.74887*10^-7 pm20-20.9433)+(5.259*10^-11 pm17-1.14698*10^-8)/ep+ep (1.86364*10^-8 pm19-54.7232)+1.34499*10^-9 pm18-2.75157*10^-7*)
vsht commented 1 year ago

Still a problem with pySecDec 1.6, I guess it must be related to issue #9

spj101 commented 3 months ago

Better late than never. I looked at this integral and wanted to share my analysis. I know you have several integrals with behaviour similar to this one, so I'm unsure to what extent the following applies generally or only to this specific integral.

Changing the propagator list to: ['-m2**2 + p1**2', '-m2**2 + (-p1 + p3)**2', '2*p3*q'] so that the Feynman parameter associated with the linear propagator is eliminated we get in total 4 sectors for this integral, which can be reduced to 2 sectors due to the symmetry of the integral.

The Symmanzik polynomials of the integral are: U^eU = ( + (1)*x0*x1)**(3*eps - 3) F^eF = ( + (m1**2)*x1*x2**2 + (m1**2)*x0*x2**2 + (m2**2)*x0*x1**2 + (m2**2)*x0**2*x1)**(1 - 2*eps) Note that U is in the denominator and leads to any divergences, while F is in the numerator.

Applying sector decomposition we find that the following monomials factor:

This means the integral has a strong singular behaviour, especially for the second sector for x1 -> 0. This results in pySecDec generating very complicated subtraction terms and also using integration-by-parts to reduce the singular behaviour in x1. This explains why the integral is slow in pySecDec especially when contour deformation is enabled and all of these steps lead to a huge expression swell.

The strange thing (that I don't yet understand) is that the integral is finite, in fact, it even starts at O(eps). Looking at the pySecDec numerical output:

    +eps^-3*(+0.0000000000000000e+00+0.0000000000000000e+00j)
    +eps^-3*(+0.0000000000000000e+00+0.0000000000000000e+00j)*plusminus
    +eps^-2*(+0.0000000000000000e+00+0.0000000000000000e+00j)
    +eps^-2*(+0.0000000000000000e+00+0.0000000000000000e+00j)*plusminus
    +eps^-1*(+7.8900801282179708e-18+0.0000000000000000e+00j)
    +eps^-1*(+5.3233755605021279e-18+0.0000000000000000e+00j)*plusminus
    +eps^0*(-1.3912426766182762e-11+0.0000000000000000e+00j)
    +eps^0*(+9.9746859433479907e-12+0.0000000000000000e+00j)*plusminus
    +eps^1*(-5.4723217606823766e+01+0.0000000000000000e+00j)
    +eps^1*(+3.0430762085916759e-10+0.0000000000000000e+00j)*plusminus
    +eps^2*(-2.0943259007907216e+01+0.0000000000000000e+00j)
    +eps^2*(+6.0967443833091971e-09+0.0000000000000000e+00j)*plusminus

We can see that the spurious 1/eps^3 and 1/eps^2 poles are analytically 0 (but this is not detected analytically by pySecDec and results in us sampling an integrand that is 0), while the 1/eps pole and eps^0 term are numerically 0, meaning we generate some non-zero integrand that nevertheless integrates to 0.

Adding a small mass to the linear propagator ['-m2**2 + p1**2', '-m2**2 + (-p1 + p3)**2', '2*p3*q-m3**2'] so that U^eU = ( + (1)*x0*x1)**(3*eps - 3) F^eF = ( + (m1**2)*x1*x2**2 + (m1**2)*x0*x2**2 + (m3**2)*x0*x1*x2 + (m2**2)*x0*x1**2 + (m2**2)*x0**2*x1)**(1 - 2*eps) and setting m3=0.0001 we actually see the same pole structure, but only the 1/eps^3 is spurious:

    +eps^-3*(+0.0000000000000000e+00+0.0000000000000000e+00j)
    +eps^-3*(+0.0000000000000000e+00+0.0000000000000000e+00j)*plusminus
    +eps^-2*(+1.3253787282080312e-09+0.0000000000000000e+00j)
    +eps^-2*(+6.4461526002390352e-25+0.0000000000000000e+00j)*plusminus
    +eps^-1*(+7.8801601152747849e-10+0.0000000000000000e+00j)
    +eps^-1*(+9.8764533617268629e-18+0.0000000000000000e+00j)*plusminus
    +eps^0*(+5.0512412030911946e-09+0.0000000000000000e+00j)
    +eps^0*(+9.9725401178695067e-12+0.0000000000000000e+00j)*plusminus
    +eps^1*(-5.4723217601432808e+01+0.0000000000000000e+00j)
    +eps^1*(+3.0431744501083905e-10+0.0000000000000000e+00j)*plusminus
    +eps^2*(-2.0943258990752952e+01+0.0000000000000000e+00j)
    +eps^2*(+6.0967355878911959e-09+0.0000000000000000e+00j)*plusminus

Next steps: