Precision loss when calculating a higher pole.

[vsht]

Hi,

I'm trying to obtain numerical results for the following integral

#!/usr/bin/env python3
if __name__ == "__main__":
    from pySecDec import sum_package, loop_package, loop_regions, LoopIntegralFromPropagators
    import pySecDec as psd

    li =     LoopIntegralFromPropagators(
    ['k2**2', 'k1**2', '-(k1*nb)', 'k2**2 + 2*gkin*k2*meta*n*u0b - k2*meta*nb*u0b - 2*gkin*meta**2*u0b**2', 'k1**2 + k3*(2*k1 >
    loop_momenta = ['k1', 'k2', 'k3'],
    powerlist = [1, 1, 1, 1, 1, 1],
    dimensionality = '4 - 2*eps',
    Feynman_parameters = 'x',
    replacement_rules = [('n**2','0'), ('nb**2','0'), ('n*nb','2'), ('gkin','1'), ('meta','1')],
    regulators = ['eps']
    )

    import os,shutil
    if os.path.isdir('loopint'):
        if input('The directory loopint already exists, do you want to overwrite it (y/n)? ') in ['y','Y','j']:
            shutil.rmtree('loopint')
        else:
            exit(1)

    loop_package(
    name = 'loopint',
    loop_integral = li,
    requested_orders = [-3],
    real_parameters = ['u0b'],
    contour_deformation = False,
    additional_prefactor = '(1)*exp(3*EulerGamma*eps)',
    decomposition_method = 'geometric'
    )

where I suspect (not 100% sure, though) that the 1/ep^3 pole should vanish.

Trying to evaluate this using Qmc

#!/usr/bin/env python3
from pySecDec.integral_interface import IntegralLibrary, series_to_mathematica, series_to_maple, series_to_sympy
import sympy as sp

li = IntegralLibrary('loopint/loopint_pylink.so')

li.use_Qmc(
)

num_params_real = [13.]
num_params_complex = []

integral_without_prefactor, prefactor, integral_with_prefactor = li(verbose = True,
real_parameters = num_params_real,
complex_parameters = num_params_complex
)

result, error = map(sp.sympify, series_to_sympy(integral_with_prefactor))

num_params_real_str = '_'.join(str(val) for val in num_params_real)
if num_params_real_str!='':
    num_params_real_str = '_' + num_params_real_str

num_params_complex_str = '_'.join(str(val) for val in num_params_complex)
if num_params_complex_str!='':
    num_params_complex_str = '_' + num_params_complex_str

print('Numerical result')
res_file = open(''.join(['numres',num_params_real_str,num_params_complex_str,'_psd.txt']),'w')
for power in [-6, -5, -4, -3]:
    val = complex(result.coeff('eps', power))
    err = complex(error.coeff('eps', power))
    out=(f'eps^{power:<2} {val: .16f} +/- {err: .16e}')
    res_file.write(out)
    res_file.write('\n')
    print(out)
res_file.close()

res_mma= series_to_mathematica(integral_with_prefactor)

res_file = open(''.join(['numres',num_params_real_str,num_params_complex_str,'_mma.m']),'w')
res_file.write(''.join(['{',res_mma[0],',']))
res_file.write(''.join([res_mma[1],'}']))
res_file.close()

res_maple= series_to_maple(integral_with_prefactor)
res_file = open(''.join(['numres',num_params_real_str,num_params_complex_str,'_maple.mpl']),'w')
res_file.write(''.join(['[',res_maple[0],',']))
res_file.write(''.join([res_maple[1],']']))
res_file.close()

I end up having huge precision issues

out.txt

Is there something I could tweak here to get a usable result?

[spj101]

Hi @vsht,

Thanks for reporting this. It seems we are really struggling with linear propagators and that they are leading to many spurious poles that complicate the subtraction.

In your case there seems to be a zero 1/eps^4 pole that is only detected upon integration (even though the integrand appears to be vanishing) and maybe a zero 1/eps^3 pole but this is less certain with the numerical uncertainty.

In your report it seems the list of propagators is truncated so we can't exactly reproduce this. Please could you send a list of all of the propagators?

[vsht]

Hi @spj101 ,

many thanks for your reply and sorry for the truncated propagator list.

I was trying to evaluate the same integral using FIESTA and initially it failed even to get the 1/ep^4 pole:

https://gitlab.com/feynmanintegrals/fiesta/-/issues/22

However, Sasha gave me a very useful piece of advice to introduce an extra regulator 'la' by hand and it
seems that with this regulator the evaluation works correctly, while the final result does not depend on the regulator.

I thought that pySecDec should be able to recognize those cases automatically, but in this particular example this
doesn't seem to be true. Anyway, if I introduce the same regulator in pySecDec, I also get that
poles up to 1/ep^2 vanish, in agreement with FIESTA.

However, FIESTA also reports that even the ep^0 term vanishes, which might be true, but is something I'd like to have checked.

Yet for some reason it seems that pySecDec doesn't want to do the higher expansions for me.

This

#!/usr/bin/env python3
if __name__ == "__main__":
    from pySecDec import sum_package, loop_package, loop_regions, LoopIntegralFromPropagators
    import pySecDec as psd

    li =     LoopIntegralFromPropagators(
    ['k2**2', 'k1**2', '-(k1*nb)', 'k2**2 + 2*gkin*k2*meta*n*u0b - k2*meta*nb*u0b - 2*gkin*meta**2*u0b**2', 'k1**2 + k3*(2*k1 + k3) - 2*gkin*(k1 + k3)*meta*n - 2*gkin*meta**2*u0b + (k1 + k3)*meta*nb*u0b', 'k2**2 + k3*(-2*k2 + k3) + 2*gkin*k2*meta*n*u0b - 2*gkin*k3*meta*n*u0b + (-k2 + k3)*meta*nb*u0b - 2*gkin*meta**2*u0b**2'],
    loop_momenta = ['k1', 'k2', 'k3'],
    powerlist = ['1 + la', '1 - la', '1 + 3*la', '1 - 3*la', '1', '1'],
    dimensionality = '4 - 2*eps',
    Feynman_parameters = 'x',
    replacement_rules = [('n**2','0'), ('nb**2','0'), ('n*nb','2'), ('gkin','1'), ('meta','1')],
    regulators = ['eps', 'la']
    )

    import os,shutil
    if os.path.isdir('loopint'):
        if input('The directory loopint already exists, do you want to overwrite it (y/n)? ') in ['y','Y','j']:
            shutil.rmtree('loopint')
        else:
            exit(1)

    loop_package(
    name = 'loopint',
    loop_integral = li,
    requested_orders = [-2,-1],
    real_parameters = ['u0b'],
    contour_deformation = True,
    additional_prefactor = '(1)*exp(3*EulerGamma*eps)',
    decomposition_method = 'geometric'
    )

works like a charm and returns

amplitude0 = + ( + ((0,0) +/- (0,0))*la^-1 + O(la^0))*eps^-2 + O(eps^-3)

in agreement with FIESTA.

However, once I try to use requested_orders = [-1,-1] I get

The directory loopint already exists, do you want to overwrite it (y/n)? y
running "sum_package" for loopint
running "make_package" for "loopint_integral"
computing Jacobian determinant for primary sector 0
total number sectors before symmetry finding: 24
total number sectors after symmetry finding (iterative): 24
total number sectors after symmetry finding (light Pak): 24
total number sectors after symmetry finding (full Pak): 24
writing FORM files for sector 1
writing FORM files for sector 2
writing FORM files for sector 3
writing FORM files for sector 4
writing FORM files for sector 5
writing FORM files for sector 6
writing FORM files for sector 7
writing FORM files for sector 8
writing FORM files for sector 9
writing FORM files for sector 10
writing FORM files for sector 11
writing FORM files for sector 12
writing FORM files for sector 13
writing FORM files for sector 14
writing FORM files for sector 15
writing FORM files for sector 16
writing FORM files for sector 17
writing FORM files for sector 18
writing FORM files for sector 19
writing FORM files for sector 20
writing FORM files for sector 21
writing FORM files for sector 22
writing FORM files for sector 23
writing FORM files for sector 24
expanding the prefactor exp(3*EulerGamma*eps)*gamma(3*eps)/(gamma(1 - 3*la)*gamma(1 - la)*gamma(la + 1)*gamma(3*la + 1)) (regulators: [eps, la] , orders: [0 2] )
 + ( + (1/3) + (0)*la + (-5*pi**2/9)*la**2)*eps**-1 + ( + (0) + (0)*la + (0)*la**2)
Traceback (most recent call last):
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polyutils.py", line 211, in _parallel_dict_from_expr_if_gens
    monom[indices[base]] = exp
          ~~~~~~~^^^^^^
KeyError: 1/eps

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/misc.py", line 505, in lowest_order
    lowest_expansion_term = sp.poly(lowest_expansion_term, variable)
                            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polytools.py", line 7299, in poly
    return _poly(expr, opt)
           ^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polytools.py", line 7272, in _poly
    result = Poly._from_expr(expr, opt)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polytools.py", line 310, in _from_expr
    rep, opt = _dict_from_expr(rep, opt)
               ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polyutils.py", line 368, in _dict_from_expr
    rep, gens = _dict_from_expr_if_gens(expr, opt)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polyutils.py", line 307, in _dict_from_expr_if_gens
    (poly,), gens = _parallel_dict_from_expr_if_gens((expr,), opt)
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polyutils.py", line 216, in _parallel_dict_from_expr_if_gens
    raise PolynomialError("%s contains an element of "
sympy.polys.polyerrors.PolynomialError: 1/eps contains an element of the set of generators.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polyutils.py", line 211, in _parallel_dict_from_expr_if_gens
    monom[indices[base]] = exp
          ~~~~~~~^^^^^^
KeyError: 1/(-5*pi**2*la**2/(9*eps) + 1/(3*eps))

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/media/Data/Projects/VS/LoopScalla/Projects/BToEtaC/Diagrams/Output/QbQubarToWQQubar/BToEtaC/3/MasterIntegrals/pySecDecFinalMasters/topology1016X01pla1mla1p3la001m3la11000/generate_int.py", line 23, in <module>
    loop_package(
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/loop_integral/loop_package.py", line 358, in loop_package
    make_package_return_value = package_generator(**LoopPackage(
                                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/make_package.py", line 375, in make_package
    sum_package(
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/code_writer/sum_package.py", line 480, in sum_package
    template_replacements = [_generate_one_term(
                             ^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/code_writer/sum_package.py", line 171, in _generate_one_term
    mp = make_package(**package_generator._asdict())
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/code_writer/make_package.py", line 2158, in make_package
    highest_prefactor_pole_orders = -np.array([lowest_order(expanded_prefactor, regulator) for regulator in regulators])
                                               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/misc.py", line 509, in lowest_order
    highest_pole = sp.poly(lowest_expansion_term**-1, variable)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polytools.py", line 7299, in poly
    return _poly(expr, opt)
           ^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polytools.py", line 7272, in _poly
    result = Poly._from_expr(expr, opt)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polytools.py", line 310, in _from_expr
    rep, opt = _dict_from_expr(rep, opt)
               ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polyutils.py", line 368, in _dict_from_expr
    rep, gens = _dict_from_expr_if_gens(expr, opt)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polyutils.py", line 307, in _dict_from_expr_if_gens
    (poly,), gens = _parallel_dict_from_expr_if_gens((expr,), opt)
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/sympy/polys/polyutils.py", line 216, in _parallel_dict_from_expr_if_gens
    raise PolynomialError("%s contains an element of "
sympy.polys.polyerrors.PolynomialError: 1/(-5*pi**2*la**2/(9*eps) + 1/(3*eps)) contains an element of the set of generators.

which looks like a bug in the code.

[vsht]

Should I open a separate bug report for the issue above? It seems unrelated to the original problem of not getting enough
precision and looks more like pySecDec for some reason forgot to expand the given term in la.