Using pySecDec only to check for the poles

[vsht]

Hi,

I have a calculation where I only need to extract the highest pole of the amplitude. This means that
a lot of master integrals get thrown away since their leading pole starts say at 1/ep or 1/ep^3 instead
of 1/ep^6.

To weed out the irrelevant integrals I'm currently using pySecDec as a mean to check which poles does
the given integral has. At this point the only information I need is that there is a pole at 1/ep^x. It's numerical
value is irrelevant for me.

However, I noticed that some integrals are so complicated, that even this kind of simple calculation takes
days and generates multiple GBs of C++-libraries. For example

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

    li =     LoopIntegralFromPropagators(
    ['k3**2', 'k2**2', 'k1**2', '(k1 + k2)**2', 'k3**2 - k3*meta*nb*u0b', 'k2**2 + 2*gkin*k2*meta*n*u0b', 'k1**2 - 2*gkin*k1*meta*n - 2*gkin*meta**2*u0b + k1*meta*nb*u0b', '(k1 + k2)**2 + k3*(-2*k1 - 2*k2 + k3) + (k1 + k2 - k3)*meta*nb*u0b', 'k1**2 + k3*(-2*k1 + k3) - 2*gkin*k1*meta*n*u0b + 2*gkin*k3*meta*n*u0b + (k1 - k3)*meta*nb*u0b - 2*gkin*meta**2*u0b**2'],
    loop_momenta = ['k1', 'k2', 'k3'],
    powerlist = [1, 1, 1, 1, 1, 1, 1, 1, 1],
    dimensionality = '4 - 2*eps',
    Feynman_parameters = 'x',
    replacement_rules = [('n**2','0'), ('nb**2','0'), ('n*nb','2')],
    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 = [-1],
    real_parameters = ['u0b', 'meta', 'gkin'],
    contour_deformation = True,
    additional_prefactor = '(1)*exp(3*EulerGamma*eps)',
    decomposition_method = 'geometric'
    )

and

#!/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 = [11., 1., 1.]
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, -2, -1]:
    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()

Since pySecDec actually needs to resolve all the poles when preparing the sector decomposition, I'm wondering if there's a better way to acheive what I want rather than just using requested_orders = [-1]. What I need is just something like

• pole at 1/ep^6 - yes
• pole at 1/ep^5 - yes
• pole at 1/ep^4 - yes
• pole at 1/ep^3 - no
• pole at 1/ep^2 - no
• pole at 1/ep^1 - no

for each integral, not more than that.

Any help would be much appreciated!

Cheers,
Vlad

[vsht]

Thanks to @magv I now know that the poles can be found inside the loopint_integral.json file, so I just need to run the generation stage.

However, I'm now struggling with some integrals where even that part doesn't really run through. For example,

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

    li =     LoopIntegralFromPropagators(
    ['k3**2', 'k2**2', 'k1**2', '(k1 + k2)**2', '-(k1*nb)', 'k3**2 - k3*meta*nb*u0b', 'k2**2 + 2*gkin*k2*meta*n*u0b', 'k1**2 - 2*gkin*k1*meta*n - 2*gkin*meta**2*u0b + k1*meta*nb*u0b', '(k1 + k2)**2 + k3*(-2*k1 - 2*k2 + k3) + (k1 + k2 - k3)*meta*nb*u0b', 'k1**2 + k3*(-2*k1 + k3) - 2*gkin*k1*meta*n*u0b + 2*gkin*k3*meta*n*u0b + (k1 - k3)*meta*nb*u0b - 2*gkin*meta**2*u0b**2'],
    loop_momenta = ['k1', 'k2', 'k3'],
    powerlist = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    dimensionality = '4 - 2*eps',
    Feynman_parameters = 'x',
    replacement_rules = [('n**2','0'), ('nb**2','0'), ('n*nb','2')],
    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 = [-1],
    real_parameters = ['u0b', 'meta', 'gkin'],
    contour_deformation = True,
    additional_prefactor = '(1)*exp(3*EulerGamma*eps)',
    decomposition_method = 'geometric'
    )

Here I never even get to the stage of having the integral library source files ready. And the memory consumption easily exceeds 1TB (!) without much happening.

Is there perhaps a better way to deal with this integral using pySecDec?

[magv]

Hi, so to document: for your example the time is spend calculating the contour deformation jacobian, which you don't need if only the poles are of interest -- disabling contour deformation should speed it up. A better solution is to use one of the sector decomposition functions directly to figure out the poles. We'll try to add a simple example for how to do this at some point.

[vsht]

Hi Vitaly,

many thanks for your help. Upon deactivating contour deformation, the above example went through really fast while the memory consumption remained pretty low.

However, there seem to be some other cases where this trick doesn't really help. For example,

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

    li =     LoopIntegralFromPropagators(
    ['k3**2', 'k1**2', 'k2*nb', '-2*gkin*meta - k1*nb', '(k1 + k2)**2 + 2*gkin*(k1 + k2)*meta*n', 'k1**2 + 2*gkin*k1*meta*n - 2*gkin*k1*meta*n*u0b', 'k2**2 + k3*(-2*k2 + k3) + (k2 - k3)*meta*nb*u0b', 'k1*(k1 - 2*k3) + k3**2 + 2*gkin*k1*meta*n - 2*gkin*k3*meta*n'],
    loop_momenta = ['k1', 'k2', 'k3'],
    powerlist = [1, 1, 1, 2, 1, 1, 4, 1],
    dimensionality = '4 - 2*eps',
    Feynman_parameters = 'x',
    replacement_rules = [('n**2','0'), ('nb**2','0'), ('n*nb','2')],
    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 = [-1],
    real_parameters = ['u0b', 'meta', 'gkin'],
    contour_deformation = False,
    additional_prefactor = '(1)*exp(3*EulerGamma*eps)',
    decomposition_method = 'geometric'
    )

Here everything goes well until the point where the log outputs writing FORM files for sector 81. From there on the memory consumption increases drastically, quickly exceeding the 32 GB RAM on my laptop.

I'm wondering if there is perhaps some trick, which could help here as well?

[magv]

Not quite sure what's going on (maybe complicated calculation of subtraction terms), but to skip writing integrands you'll need to use sector decomposition functions directly.

[vsht]

An example would be really nice. In the meantime it seems that FIESTA can handle my problematic integrals more efficiently, even evaluating them numerically. Since there are only O(20) of them, I hope that I can avoid diving into the pySecDec internals for now ;)

Even though I usually prefer to use pySecDec, 2TB RAM usage on a preparation stage with contour deformation turned off looks really weird to me.

[j-schlenk]

Hi,
here is a function that returns the maximal pole order for a given integral (here the one-loop box) without writing out integrands:

#!/usr/bin/env python3
import pySecDec as psd
import numpy as np

if __name__ == "__main__":

    def pole_order(li, includegamma=True):
        vertices = np.concatenate((li.U.expolist, li.F.expolist))
        polytope = psd.polytope.Polytope(vertices = vertices)
        polytope.complete_representation()
        facets = polytope.facets
        dimension = li.dimensionality.subs([(regulator,0) for regulator in li.regulators])
        intpoint = li.powerlist + [dimension/2]

        incidence_lists = polytope.vertex_incidence_lists()
        transformation = polytope.facets.T[:-1]
        poleorder = facets.dot(intpoint)
        poleorderdict = {tuple(fac): i for i,fac in zip(poleorder,(transformation.T).tolist())}

        poleorderglobal = 0
        for vertex, cone_indices in zip(incidence_lists.keys(),incidence_lists.values()):
            cone = transformation[:,cone_indices].T
            if len(cone_indices) != len(li.powerlist):
                triangular_cones = psd.polytope.triangulate(cone)
            else:
                triangular_cones = np.array([cone])
            poleordersector = max(sum( 1 for x in [poleorderdict[tuple(facet)] for facet in cone] if x <=0) for cone in triangular_cones.tolist())
            poleorderglobal = max([poleordersector, poleorderglobal])
            
        # pole from Gamma[N-L*d/2] prefactor in Feynman parametrisation
        if sum(intpoint[:-1])-intpoint[-1]*(li.L)<=0 and includegamma == False:
            poleorderglobal = poleorderglobal - 1
        #remove spurious expoU<0 pole from 1/Gamma[-N+(L+1)d/2] in Lee-Pomeransky representation
        if -sum(intpoint[:-1])+intpoint[-1]*(li.L+1)<=0:
            poleorderglobal = poleorderglobal - 1

        return poleorderglobal
  
    # Define the integral
    li = psd.loop_integral.LoopIntegralFromPropagators(
        loop_momenta = ['k1'],
        external_momenta = ['p1','p2','p3'],
        propagators = ['k1**2','(k1+p1)**2','(k1+p1+p2)**2','(k1+p1+p2+p3)**2'],
        powerlist = [1,1,1,1],
        dimensionality = '4-2*eps',
        replacement_rules = [
                        ('p1*p1', 0),
                        ('p2*p2', 0),
                        ('p3*p3', 0),
                        ('p1*p2', 's/2'),
                        ('p2*p3', 't/2'),
                        ('p1*p3', '-s/2-t/2')
                    ]
    )

    print(pole_order(li,includegamma=True))

The function should reproduce the pole order that is found by pySecDec.
This means that the returned pole order can be too high if there are spurious poles (for example the integral in issue #27).
The flag includegamma determines if the pole from the Gamma function $\Gamma(N-Ld/2)$ in the prefactor of the Feynman parametrisation is included in the pole order counting.

[spj101]

In release 1.6.5 we have now introduced the function loop_integral.leading_order().
Hopefully this returns the information you are looking for.