Some questions regarding the usage of HyperInt

[Vladyslav Shtabovenko]

Hello,

I apologize in advance if it is not the right channel to ask these kind of questions, but I’m wondering how one should approach the issue when an integration in a Feynman parameter breaks down with the error message

Divergence at x[2] = infinity of type ln(x[2])

As most other users, I’m employing HyperInt to calculate loop integrals (in this case 3-loop propagator-type with several massive lines) via direct Feynman parametrization. Although in many cases things simply work out of the box (which is pretty cool), there are also situations (like here), where I get stuck on the error messages that are not entirely clear to me. In such cases I try to check with the manual of play with my Feynman parametrization (e.g. by joining the propagators not all at once, but sequentially).

Long story short, I’m attaching an example where the error pops up. It would be enough for me to know, whether the issue can be resolved in a trivial way (i.e. by changing some of the HyperInt parameters) or whether it is an artefact of a bad choice of Feynman parametrization.

‌

‌

‌

‌

[Erik Panzer]

Dear Vladyslav,

thanks for reporting your problem! I had a look, this is what’s happening: With

eval(_hyper_divergences);

or, in this particular case (as indicated by the error message),

C:=_hyper_divergences[x[2]=infinity,ln(x[2])];

you can look at the coefficient of ln(x[2]), causing a suspected divergence at x[2]->infinity. You will see logarithms and dilogarithms with 3rd roots of unity in the arguments. Since HyperInt’s period tables only come with multiple zeta values and alternating sums, HyperInt cannot simplify this transcendental number. It is in fact zero, which you can see with

evalf(C);
evalf[50](C);

which gives zero up to the given precision. So there is in fact no divergence, but HyperInt aborts, because it cannot symbolically establish that this number C is zero.

The quick solution is to simply put

_hyper_abort_on_divergence:=false;

which tells HyperInt to continue. This is precisely made for such situations. It will still compute the potentially divergent terms, that it cannot identify with zero, and store them in _hyper_divergences. This allows you to check afterwards if there really was a divergence or not.

Alternatively, if you have a lot of integrals with 3rd roots of unity to process, it may be better to create a table that reduces the polylogarithms of low weight with such roots as arguments, down to a basis. (e.g. using the table from https://arxiv.org/abs/1512.08389).

I hope this helps - please let me know if that resolves your problem.

With best wishes,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

many thanks for your detailed answer! This is really helpful to understand the logic behind HyperInt’s error messages.
In principle, it would be already sufficient for me to have an analytic results in terms of functions that can be evaluated using
HyperInt or PolyLogTools. The reduction to a basis is, of course, always nice, but is probably not necessary for this project.

May I also ask, how should one handle _hyper_splitting_field when extending the field beyond simple cases such as I or a single RootOf?
For example, I noticed that _hyper_splitting_field := {I,RootOf(_Z^2 + _Z + 1)}; apparently does not represent a valid input, since it leads to

Error, (in factors) 2nd argument, {I, RootOf(_Z^2+_Z+1)}, is not a valid algebraic extension

Furthermore, it is not entirely clear to me how to make HyperInt consider also x-variables when factoring terms. In the following small example

_hyper_splitting_field := {};_hyper_abort_on_divergence:=false;_hyper_algebraic_roots:=true;forgetAll():

expr:=[[12*(x[3]^2x[6]^2+2x[3]^2x[6]+2x[3]x[6]^2+x[3]^2+2x[3]x[6]+x[6])/(x[6]+1)/(x[3]x[6]+x[3]+x[6])^2/x[6], [[-(x[3]^2x[6]+x[3]^2+2x[3]*x[6]+x[3]+x[6])/(x[3]*x[6]+x[3]+1), 0, -x[6]/(x[6]+1)]]]];

hyperInt(expr,x[6])

hyperInt apparently doesn’t like x[3]^2_x[6]+x[3]^2+2_x[3]*x[6]+x[3]+x[6] : replacing this numerator in the second slot by unity as in

expr:=[[12*(x[3]^2x[6]^2+2x[3]^2x[6]+2x[3]x[6]^2+x[3]^2+2x[3]x[6]+x[6])/(x[6]+1)/(x[3]x[6]+x[3]+x[6])^2/x[6], [[-(1)/(x[3]*x[6]+x[3]+1), 0, -x[6]/(x[6]+1)]]]];

makes everything work. However, it seems that writing

_hyper_splitting_field := {RootOf(x[3]^2x[6]+x[3]^2+2x[3]*x[6]+x[3]+x[6],x[6])};

is not a correct way to solve this issue. In fact, in the case of x[3]^2_x[6]+x[3]^2+2_x[3]*x[6]+x[3] (i.e. with a single x[6] removed) the integration also succeeds and contains explicit Root(x[3]^2_x[6]^2+2_x[3]^2*x[6]+x[3]_x[6]^2+x[3]^2+2_x[3]*x[6]+x[3]-x[6], x[6]). In this sense, I’m a bit puzzled with the current behavior.

However, in this special case it appears that even if I can somehow get past the integration in x[6], the last integration in x[3] still would not work

expr:=[[12*(x[3]^2x[6]^2+2x[3]^2x[6]+2x[3]x[6]^2+x[3]^2+2x[3]x[6]+x[6])/(x[6]+1)/(x[3]x[6]+x[3]+x[6])^2/x[6], [[-(x[3]^2x[6]+x[3]^2+2x[3]*x[6]+x[3])/(x[3]*x[6]+x[3]+1), 0, -x[6]/(x[6]+1)]]]];
ex2:=hyperInt(expr,x[6]);
hyperInt(ex2,x[3])

Perhaps I should really try to find a better parametrization for this integral.

‌

[Erik Panzer]

Dear Vladyslav,

splitting fields should be a list or set of radicals, or a list or set of RootOf’s, see https://www.maplesoft.com/support/help/maple/view.aspx?path=factor

In particular, you are not supposed to mix RootOf’s and other radicals. So your

_hyper_splitting_field := {I,RootOf(_Z^2 + _Z + 1)};

is not allowed, but instead you could write

_hyper_splitting_field := {I,sqrt(3)};

because the third roots of unity are rational linear combinations of I and sqrt(3). If you prefer RootOf’s, you could write the same field as

_hyper_splitting_field := {RootOf(_Z^2+1),RootOf(_Z^2 + _Z + 1)};

or

_hyper_splitting_field := {RootOf(_Z^2+1),RootOf(_Z^2-3)};

Regarding your other question, the problem is that HyperInt only works for rational functions, and accounts for RootOf’s only within fibrationBasis. The problem with your integral is that it is not dlog, because the rational prefactor of expr has a term

(24*x[3]+12)/(x[3]*x[6]+x[3]+x[6])^2

in its partial fraction decomposition wrt x[6]. The integration of x[6] therefore produces a term, which still needs to be integrated, that has prefactor (24*x[3]+12)/(x[3]*x[6]+x[3]+x[6])/(1+x[3]), but which multiplies a derivative of a hyperlogarithm with RootOf’s. This brings the RootOf’s in the prefactor, and HyperInt just does not deal with that.

‌

One could consider a workaround and extending HyperInt to accept arbitrary RootOf’s, however, this is not trivial and can be very slow, and moreover the main point is that this would already happen at your penultimate integration. The final integration would then be impossible, because the integrand would not have rational, but RootOf(…x[3]…)-dependent arguments in the hyperlogarithms.

‌

So it seems that you try to integrate a non-linearly reducible integral. Have you checked the polynomial reduction?

Best,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

I see, many thanks for the clarification. I have to admit that I’m usually using Mathematica and the only reason I started to look at Maple is
HyperInt. So I’m still 100% comfortable with the Maple syntax, but I’m doing my best to improve.

You are of course right regarding the linear reducibility, I should have payed more attention to that from the very beginning . In my case (the original attached notebook), using checkIntegrationOrder on the 4 of the 6 integration variables leaves me with

{x[2] + x[3], 2*x[3] + x[2], x[2]^2 + x[2]*x[3] + x[3]^2}

which clearly asks for a variable transformation. I managed to come up with xp[2]->x[2]+x[3], xp[3]->x[2]-x[3] which leads to

{x[2] - x[3], x[2] + x[3], x[2] + x[3]/3, x[2]^2 + x[3]^2/3}

This way I can integrate O(ep^0) and O(ep^1) pieces directly. The latter requries _hyper_splitting_field := {RootOf(_Z^2 + 1/3)}: and the result seemingly cannot be expressed in terms of MZVs alone, but this is already something. The issue with (x[3]*x[6]+x[3]+x[6]) still hits me at O(ep^1), e.g.

_hyper_algebraic_roots := true:
_hyper_splitting_field := {RootOf(_Z^2 + _Z + 1)}:
_hyper_ignore_nonlinear_polynomials:=false;
_hyper_abort_on_divergence:=true;
_hyper_verbosity:=2;
forgetAll():
hyperInt([[-96*(x[3]*x[6] - x[3] + x[6] + 1)/((x[3]^2 - 2*x[3]*x[6] - 2*x[6] - 1)*(x[3]^2 - 4*x[6] - 1)*(x[3] - 2*x[6] - 1)), [[-(x[3]^2 - 4*x[6] - 1)/(2*x[3] - 4*x[6] - 2), -(x[3]^2*x[6] - x[3]^2 + 3*x[6] + 1)/(2*(x[3]*x[6] - x[3] + x[6] + 1)), x[6]*(x[3] - 1)/(x[3]*x[6] - x[3] + x[6] + 1)]]]],x[6]);

But I suppose that the general strategy to deal with it is to introduce another transformation that simplifies denominators in x[3] and x[6] and so on.

By the way, I’m very grateful for your explanations. HyperInt is an enormously powerful tool, but it is not always dimplr to learn how to use it in the most efficient way.

Cheers,
Vladyslav

‌

[Erik Panzer]

Dear Vladyslav,

I had a closer look at your file now. Actually, since you are projective, there are only 5 integrals to do, so your integral is linearly reducible! No need for any reparametrizations. I do not understand where that integral in your last post comes from, but I attach a file that computes what you might be after.

read HyperIntFile:
psi:=x[1]*x[2]*x[3] + x[1]*x[2]*x[4] + x[2]*x[3]*x[4] + x[2]*x[3]*x[5] + x[2]*x[4]*x[5] + x[1]*x[2]*x[6] + x[1]*x[3]*x[6] + x[2]*x[3]*x[6] + x[1]*x[4]*x[6] + x[3]*x[4]*x[6] + x[2]*x[5]*x[6] + x[3]*x[5]*x[6] + x[4]*x[5]*x[6];
phi:=x[1]*x[2]^2*x[3] + x[1]*x[2]*x[3]^2 + x[1]*x[2]^2*x[4] + x[1]*x[2]*x[3]*x[4] + x[2]^2*x[3]*x[4] + x[2]*x[3]^2*x[4] - x[1]*x[2]*x[3]*x[5] + x[2]^2*x[3]*x[5] + x[2]*x[3]^2*x[5] - x[1]*x[2]*x[4]*x[5] + x[2]^2*x[4]*x[5] + x[1]*x[2]^2*x[6] + x[1]*x[2]*x[3]*x[6] + x[2]^2*x[3]*x[6] + x[1]*x[3]^2*x[6] + x[2]*x[3]^2*x[6] + x[1]*x[3]*x[4]*x[6] + x[3]^2*x[4]*x[6] - x[1]*x[2]*x[5]*x[6] + x[2]^2*x[5]*x[6] - x[1]*x[3]*x[5]*x[6] + 2*x[2]*x[3]*x[5]*x[6] + x[3]^2*x[5]*x[6] - x[1]*x[4]*x[5]*x[6];
F:=((x[4]*x[5] + x[3]*(x[4] + x[5]))*x[6] + x[1]*((x[3] + x[4])*x[6] + x[2]*(x[3] + x[4] + x[6])) + x[2]*(x[5]*(x[4] + x[6]) + x[3]*(x[4] + x[5] + x[6])))^(-2 + 4*ep)*(x[3]^2*(x[4] + x[5])*x[6] + x[1]*((x[3] + x[4])*(x[3] - x[5])*x[6] + x[2]^2*(x[3] + x[4] + x[6]) + x[2]*(x[3] - x[5])*(x[3] + x[4] + x[6])) + x[2]*x[3]*(2*x[5]*x[6] + x[3]*(x[4] + x[5] + x[6]))+x[2]^2*(x[5]*(x[4] + x[6]) + x[3]*(x[4] + x[5] + x[6])))^(-3*ep);
findDivergences(F);
F2 := factor(dimregPartial(F, {x[2],x[3],x[4],x[6]}, 2*ep)): findDivergences(F2);
F3 := factor(dimregPartial(F2, {x[2],x[6]}, ep)): findDivergences(F3);
L := table():
S := irreducibles({phi, psi});
L[{}] := [S, {S}]:
cgReduction(L):
suggestIntegrationOrder(L);
vars := [x[1], x[5], x[4], x[6],x[2],x[3]];
checkIntegrationOrder(L, vars[1..4]);
temp:=hyperInt(eval(coeff(series(F3, ep = 0), ep, 0), vars[-1] = 1), vars[1 .. -3]):
_hyper_splitting_field:=RootOf(_Z^2+_Z+1):
_hyper_abort_on_divergence:=false:
Result:=hyperInt(temp,vars[-2]);
simplify(fibrationBasis(Result));

Please let me know if this is what you intended.

‌

Also, if you have any ideas of which sections in the manual could need improvement, please tell me.

‌

By the way, your integral seems to be ill-defined. The final answer has a term +-15/2*i*Pi, whose sign depends on whether x[2] gets a +i*epsilon or -i*epsilon during the last integration. This means that the integrand of the last integration is discontinuous on the positive real line. Do you expect such a multi-valued/ill-defined result?

‌

With best wishes,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

sorry for the late reply. Thanks for having a look at my example. In fact, I wrote a similar code after you explained me the trick with
_hyper_abort_on_divergence:=false; Even though one can get the O(ep^0) piece of the result, the Feynman parametrization actually contains
a GAMMA(3*ep) so that one has to go at least one order in ep higher. In this case already

temp:=hyperInt(eval(coeff(series(F3, ep = 0), ep, 1), vars[-1] = 1), vars[1 .. -3]):

ends up with

Error, (in linearFactors) x[2]^2x[6]^2+x[2]^3+2x[2]^2*x[6]+x[2]x[6]^2+x[2]^2+2x[2]*x[6]+x[6]^2 does not factor linearly in x[6]

This was the reason for my other inquiry and from your reply I understood that this is not something that can be easily resolved. Hence came the idea
to try to find a better parametrization. The one that I initially suggested is not correct, perhaps I should look for a better one.

I’m not quite sure what you mean with “ill-defined “. The integral in question (one of many that I’d like to calculate) -
{l1^2, l2^2, (l1 + q)^2, (l2 + l3 + q)^2, l3^2 - m^2, (l1 - l2)^2 - m^2} is an on-shell (q^2 =m^2) 3-loop integral with 2-massive lines.
As you correctly observed, the letters will be somewhat weird, but I see no intrinsic issues with the integral itself.
The integral is of course projective, otherwise I would have projectivized it before passing to HyperInt.

Perhaps I’m a too simple-minded person, but the appearance of delta[x[2]] doesn’t necessarily make me worried.
At the end of the day, all analytic results are always cross-checked using FIESTA/secDec and the numerics can be also used to fix the signs
of the deltas.

In fact, I become more worried when after applying evalf to HyperInt’s result containing RootOf(_Z^2 + _Z + 1), I don’t get neither the deltas
nor the expected imaginary part. But of course, there is always a chance that I did something wrong inbetween. I can also provide the code, if needed
(this involves another integral).

In general, the reason for my questions is not to calculate some particular integrals, but rather to better understand how to use HyperInt in my calculations.

I think that the manual is already very good, but perhaps it could provide more examples for nonstandard case (like the one in my original post),
where it may not be completely obvious how to obtain a final result despite of the error messages. I mean, if the integral is already projective and well
behaved, then applying HyperInt is extremely simple and one can blindly copy paste commands from the very first examples in the manual. However, once
things do not work out this way, one is left wondering whether one should try harder or just switch to another method to calculate the integral.

As a practical suggestion, perhaps one could add a routine that projectivizes an integral that is not already projective. I have my own code for that,
but since this is a rather simple operation I was wondering why something like that cannot be integrated directly into HyperInt.

[Vladyslav Shtabovenko]

Dear Erik,

In fact, I become more worried when after applying evalf to HyperInt’s result containing RootOf(_Z^2 + _Z + 1), I don’t get neither the deltas
nor the expected imaginary part. But of course, there is always a chance that I did something wrong inbetween. I can also provide the code, if needed
(this involves another integral).

please disregard the above part of my message, I had a bug in the code assembling the whole expression. For the integral in question (not the one mentioned in my very first post) HyperInt produces a correct results with a `delta[x[6]]` indicating the +/- 1 uncertainty in the imaginary part. But since this can be resolved by comparing to a numerical evaluation, it is not really an issue for me.

There are actually two points that came to mind and would constitute a nice extension of the manual

• Explain how to extend fibrationBasis to support new letters (along the lines of what you suggested to me regarding arXiv:1512.08389
• Show how one can obtain an analytic result without such an extension (but when the letters are known) by using PSLQ

I guess that would be something interesting for people using HyperInt in loop calculations.

Cheers,
Vladyslav

[Erik Panzer]

Dear Vladyslav,

thank you for these suggestions! Indeed, the support and documentation for non-harmonic transcendentals is very poor at the moment. I’ll expand that in the manual (hopefully soon). For the purpose of reducing the kind of numbers that emerge from your integrals to a basis, I would suggest to consider also Oliver Schnetz’es HyperlogProcedures code. It uses Francis' Brown’s motivic decomposition algorithm, and is very efficient.

‌

Regarding the problem with your integral due to

Error, (in linearFactors) x[2]^2x[6]^2+x[2]^3+2x[2]^2*x[6]+x[2]x[6]^2+x[2]^2+2x[2]*x[6]+x[6]^2 does not factor linearly in x[6]

Note that the polynomial reduction does not depend on the order in epsilon. So, the fact that the polynomials {U,F} are linearly reducible, implies that you can integrate all orders in the epsilon expansion. Hence, any non-linear polynomials that appear must be spurious (i.e. occur for individual summands, but cancel out completely when combined to the full integrand). So you can simply put

_hyper_ignore_nonlinear_polynomials := true;
temp := hyperInt(eval(coeff(series(F3, ep = 0), ep, 1), vars[-1] = 1), vars[1 .. -3]);

to get the result. The best way to exploit the polynomial reduction is explained in the documentation and in the manual (using _hyper_restrict_singularities and _hyper_allowed_singularities), but if efficiency is not a problem, you can simply tell HyperInt to ignore any non-linear polynomials, as shown above.

‌

Does this solve your problem?

‌

With best wishes,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

first of all I’d like to express my sincere gratitude for your numerous advices. By combining all the hints you gave me
I was able to calculate the most complicated of my 3-loop integrals with HyperInt with very little effort. For most of those
integrals I also have a numerical result from pySecDec, which was handy to fix the values of delta[x[i]] to be +1 or -1.
Needless to say that the analytic results are in perfect agreement with the output of pySecDec. We will certainly acknowledge
your kind help in the corresponding paper.

Essentially, I have now only two minor questions left.

1. The results for two of my integrals (3-loops, 6 lines, still on-shell) contain a zeta[1,-3] in the real part. Do I understand it correctly that
   it signals some noncancellation of “poles” in the regularization of HyperLogarithms, or is there a simpler explanation? One of the offending
   terms appearing at O(ep^1) looks like [[-24, [[0, -1/2, 0, -1]]]] and becomes

-219/5*zeta[2]^2+63*ln(2)*zeta[3]-36*zeta[1,-3]-12*ln(2)^2*zeta[2]-3*ln(2)^4

upon applying fibrationBasis to it. The matter is not too pressing, since we actually need only the imaginary part (and that one is obviously correct), but I’d like
to better understand how to interpret such output. I can, of course, provide the full calculation if you need it.

2. If I understood you correctly, in order to express my analytic results (that are currently written as Hlogs with many RootOf and look e.g. as follows

-Hlog(1,[0, 1/(1+RootOf(_Z^2+_Z+1))])-Hlog(1,[0, -1/RootOf(_Z^2+_Z+1)])-2*Hlog(1,[0, 1/(1+RootOf(_Z^2+_Z+1)), 1])-1/6*ln(1+RootOf(_Z^2+_Z+1))^3-1/6*ln(-(-1-RootOf(_Z^2+_Z+1))/(1+RootOf(_Z^2+_Z+1)))^3-2*Hlog(1,[0, -1/RootOf(_Z^2+_Z+1), 1])-1/6*ln(-RootOf(_Z^2+_Z+1))^3-1/2*ln(1+RootOf(_Z^2+_Z+1))^2*ln(-(-1-RootOf(_Z^2+_Z+1))/(1+RootOf(_Z^2+_Z+1)))-1/2*ln(1+RootOf(_Z^2+_Z+1))*ln(-(-1-RootOf(_Z^2+_Z+1))/(1+RootOf(_Z^2+_Z+1)))^2+1/2*ln(-(-1-RootOf(_Z^2+_Z+1))/(1+RootOf(_Z^2+_Z+1)))^2+1/2*ln(1+RootOf(_Z^2+_Z+1))^2+1/2*ln(-RootOf(_Z^2+_Z+1))^2+2*Zeta(3)-11/360*Pi^4-Hlog(1,[-1/RootOf(_Z^2+_Z+1), 0, 1, 1])-Hlog(1,[0, -1/RootOf(_Z^2+_Z+1), 1, 1])-Hlog(1,[-1/RootOf(_Z^2+_Z+1), 1, 0, 1])-Hlog(1,[1/(1+RootOf(_Z^2+_Z+1)), 0, 1, 1])-Hlog(1,[0, 1/(1+RootOf(_Z^2+_Z+1)), 1, 1])-Hlog(1,[1/(1+RootOf(_Z^2+_Z+1)), 1, 0, 1])-4*Hlog(1,[1/(1+RootOf(_Z^2+_Z+1)), 0, 0, 1])+Hlog(1,[0, 0, 1/(1+RootOf(_Z^2+_Z+1)), 1])-Hlog(1,[0, 0, 0, 1/(1+RootOf(_Z^2+_Z+1))])-Hlog(1,[0, 0, 0, -1/RootOf(_Z^2+_Z+1)])-4*Hlog(1,[-1/RootOf(_Z^2+_Z+1), 0, 0, 1])+Hlog(1,[0, 0, -1/RootOf(_Z^2+_Z+1), 1])-ln(1+RootOf(_Z^2+_Z+1))-ln(-(-1-RootOf(_Z^2+_Z+1))/(1+RootOf(_Z^2+_Z+1)))+ln(1+RootOf(_Z^2+_Z+1))*ln(-(-1-RootOf(_Z^2+_Z+1))/(1+RootOf(_Z^2+_Z+1)))+I*Pi*Zeta(3)-RootOf(_Z^2+_Z+1)*ln(-(-1-RootOf(_Z^2+_Z+1))/(1+RootOf(_Z^2+_Z+1)))-RootOf(_Z^2+_Z+1)*ln(1+RootOf(_Z^2+_Z+1))+RootOf(_Z^2+_Z+1)*ln(-RootOf(_Z^2+_Z+1))+Hlog(1,[0, 0, -1/RootOf(_Z^2+_Z+1)])-4*Hlog(1,[-1/RootOf(_Z^2+_Z+1), 0, 1])-1/6*ln(-RootOf(_Z^2+_Z+1))*Pi^2-1/6*ln(1+RootOf(_Z^2+_Z+1))*Pi^2-1/6*ln(-(-1-RootOf(_Z^2+_Z+1))/(1+RootOf(_Z^2+_Z+1)))*Pi^2+2*Zeta(3)/ep+Hlog(1,[0, 0, 1/(1+RootOf(_Z^2+_Z+1))])-4*Hlog(1,[1/(1+RootOf(_Z^2+_Z+1)), 0, 1])+I*Pi+1/6*Pi^2-1/6*I*Pi^3;

) in terms of simpler expressions, I could try to create a weight reduction table for HyperInt or use HyperlogProcedures. I also thought that since HyperInt can evaluate
my expressions with arbitrary precision, PSLQ might be a good option, provided that I can figure out that complete alphabet. What would be the most promising approach in your opinion?

I would be actually quite eager to learn the procedure of creating reduction tables for new letters, but I’m not sure how complicated this is for real life problems.

PS: Nochmal vielen Dank für deine Hilfe! Ist bestimmt alles nur Kinderkram für dich, aber als jemand, der ansonsten nicht so viel mit der Berechnung der Masterintegrale
am Hut hat, finde ich deine Erläuterungen enorm hilfreich!

Cheers,
Vladyslav

[Vladyslav Shtabovenko]

Sorry, the first question makes no sense. I obviously forgot that the MZV convention used in HyperInt has the arguments list reversed w.r.t the convention
used in HPL/PolyLogTools (facepalm). Hence zeta[1,-3] is MZV[{-3,1}] ~ 0.08778 and the real part is in perfect agreement with pySecDec.

I’m currently looking into the reduction of my HyperLogs to simpler function, but any kind of hints (even sketchy ones) would be very much appreciated.

Cheers,
Vladyslav

[Erik Panzer]

Dear Vladyslav,

‌

I am glad that you managed to get through the calculation as such! Thinking about the reductions, I can write a bit more about using HyperlogProcedures in a bit.

The most convenient solution, in particular for your application that seems to involve only small weights, would be inside HyperInt via a reduction table. The manual and documentation explain how this needs to be done, there should be an example for 4th roots of unity that hopefully illustrates how this goes.

As for the actual source of the reductions, upon further reflection, your result seems to lie in the subspace of multiple Deligne values. In particular, it does not mix all 6th roots of unity, but instead only has letters 0, 1 and the primitive 6th roots of unity. David Broadhurst has provided a handy table for these up to weight 11 in https://arxiv.org/abs/1409.7204.

For example, the first Hlog in your expression is

Hlog(1,[0, 1/(1+RootOf(_Z^2+_Z+1))]) = -Z(AD)

in Broadhurst’s notation, and this encodes the dilogarithm (use convert(…,Mpl))

Hlog(1,[0, 1/(1+RootOf(_Z^2+_Z+1))]) = -polylog(2,1/2+I/2*sqrt(3))

whose imaginary part is the generator I2 in Broadhurst’s tables.

This way, you can look up all Hlog’s in your expression in Broadhurst’s tables and reduce it to an expression involving only Pi, Zeta(3), the dilogarithm above, and products thereof (but I could well have made a mistake).

If you process Broadhurst’s files into HyperInt-compatible format (see the 4th roots example), you could load them into HyperInt and they would be substituted automatically in fibrationBasis(), just like MZVs. This would be the preferred solution. Unfortunately, I have not yet worked myself with these multiple Deligne values, so I have no Maple version of the files prepared yet.

With best wishes,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

as it turned out, my integrals actually contain GPL with more letters than just A,B,C and D from the paper David Broadhurst. Using the reduction tables
from arXiv:1512.08389 and doing quite some PSLQ to get rid of particular combinations of the basis elements I at least managed to simplify the imaginary
parts of the integrals (which is what we are mainly after) to combinations of Pi, Zeta’s, Log(2), Log(3) and a Clausen(2,Pi/3). Claude Duhr’s PolyLogTools
actually implements some of the relations between GPLs with arguments that are sixth root of unity, but only up to weight 3, while I had to deal with cases
with up to weight 6.

All in all, it was a lot of fun and I feel that I really learned a lot in the past two weeks. And for sure without your help I wouldn’t have hardly made so much
progress in so little time.

One particular issue I noticed is that the dimRegPartial-method of regularizing divergences can significantly slow down the integration for integrals
that contain multiple singularities in different combinations of Feynman parameters. Now that I still have to calculate several 5-line 3-loop integrals
up to O(ep^2) this property starts to become a limitation, since the evaluations literally take many hours.

I suppose that the preferred solution is to use the methods outlined in your papers with Robert Schabinger and Andread von Manteuffel
(arXiv:1411.7392, arXiv:1510.06758). Perhaps this would also be something worth adding to the manual.

In any case, I think that all of my questions regarding the usage of HyperInt were nicely answered by you and I am very grateful for that.
I also hope that the information provided here can be useful for other users of HyperInt.

I also apologize for misusing the issue tracker as a sort of “forum”. To my knowledge, Bitbucket currently doesn’t offer anything similar to
GitHub Discussions, so that it’s sometimes difficult to find the correct communication channel. Perhaps you could also add some relevant
information with respect to that to README.md, i.e. explain how should people proceed if they questions about HyperInt that are not bug reports.

Thanks again!

Cheers,
Vladyslav

[Vladyslav Shtabovenko]

All my questions have been answered by the package author.

[Erik Panzer]

Dear Vladsylav,

thank you for your great suggestions! I added some HyperInt output regarding _hyper_abort_on_divergence and updated the documentation accordingly, which may help with one of the initial issues that you encountered.

Regarding the periods, I added references for the 6th roots of unity tables to section 4.1 in the documentation. It would be nice to turn those tables into HyperInt-ready Maple files and provide them as a download for other users with 6th roots of unity.

I also added a comment to the end of section 4, referring to the motivic decomposition algorithm. However, I still need to add an actual example, but at least it mentions the possibility and the great utility of HyperlogProcedures.

I suppose that the preferred solution is to use the methods outlined in your papers with Robert Schabinger and Andread von Manteuffel
(arXiv:1411.7392, arXiv:1510.06758). Perhaps this would also be something worth adding to the manual.

Indeed! If you have the IBP under control, switching to a finite basis would be the easiest starting point for the integration, avoiding all those dimregPartial’s. I took a note to add a corresponding reference to the documentation soon.

If you cannot deal with the IBP, sometimes it can help if you resolve the divergences in a different order - depending on the order of dimregPartial’s, the expression sizes can differ vastly.

There is also a pending improvement possible for dimregPartial, which is alas not implemented yet.

Perhaps you could also add some relevant
information with respect to that to README.md, i.e. explain how should people proceed if they questions about HyperInt that are not bug reports.

With the current setup, I think the issue tracker, in the way you used it, is a rather effective way to communicate these things. Only that after resolution, posts tend to disappear and be hard to find, so it might be helpful to keep a list with links to potentially helpful issue discussions on the wiki (of course that would be a very cheap replacement of a forum, but since the number of these discussions can be counted on one hand, it may just about do).

If you have any other ideas or suggestions, I’d be happy to hear.

‌

Best,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

thank you for all the improvements and hints and especially for adding a direct conversion to the notation used
in HyperLogProcedures. Using Convert(res,f23) I was able to simplify my integrals straightforwardly
to very compact expressions in just few minutes, although in a few cases the Convert procedure
generated an error of type “Error, (in evala/preproc3) floats not handled”.
A simplest code sample that reproduces this behavior would be e.g.

Convert(3*i[1, 1/(1 + RootOf(_Z^2 + _Z + 1)), 1, 2, 1, 1, 0],f23)

For the sake of completeness, I’m also attaching full expressions for the offending integrals (see above).
I suppose it is because the letter 2 is spurious here and doesn’t belong to f23. In fact, by splitting such GPLs
into real and imaginary parts one can show that they cancel in the final result. Still, I’m not sure if the output of
HyperLogProcedures in this case is ok, or if there is a need to contact the author.

I’m also somewhat surprised that the conversion from the f23- back to the i-basis seems to require almost a minute even for very simple
expressions such as

Convert(5*I*sqrt(3)*f23[2]/(2*ep),i)

Probably it would be easier for me to do this via applyrule since HyperLogProcedures is somehow incredibly slow here.

I’m also wondering if HyperInt could possibly handle the conversion of the i-functions back into Hlogs or Mpls? Not that one couldn’t
do it via applyrule, but this still could be quite handy.

Cheers,
Vladyslav

[Erik Panzer]

Dear Vladyslav,

I am not sure what the origin of this error is. Please contact Oliver Schnetz.

For what it’s worth, I fed your numbers into an older version (0.2) of HyperlogProcedures, and there I get results without errors:

Convert(3*i[1, 1/(1 + RootOf(_Z^2 + _Z + 1)), 1, 2, 1, 1, 0],f23)
= 5/108*I*Pi^3*f23[2]-117/64*f23[2,3]+59/384*Pi^2*f23[3]+531/128*f23[3,2]+3/2*I*Pi*f23[4]-2775/256*f23[5]-179/50544*I*Pi^5;


resFin 
= 1/1080*Pi^4*(81*I*Pi*delta[x[2]]-64)+(7/2*Pi^2-3/8*Pi^2*delta[x[2]]^2)*f23[3]-\
63/8*I*Pi*f23[4]-27/40*f23[5]-4/135*Pi^4/ep


resFin2 
= 1/8*I*Pi^3-383/6480*Pi^4+15/2*I*Pi-25/8*Pi^2+95/2-5/4*I*Pi*f23[-1,2]+21/8*f23[
-1,3]-1/4*I*Pi^3*f23[-1]+1/2*I*Pi*f23[2]+(-219/16+45/32*I*Pi)*f23[3]+(5/2*I*Pi
-73/72*Pi^2+55/6)/ep+(1/2*I*Pi+3/2)/ep^2+1/6/ep^3

If you can check that they work numerically, they might be correct. From the fact that they run in version 0.2 but not in 0.4, Oliver might be able to quickly figure out the origin of the issue.

Regarding the back-conversion from f’s to i’s, I think it is expected that this is rather slow (the whole idea is to simplify i’s by rewriting them as f’s, and to go the other way, one is essentially forced to build up more or less a full table again, if I am not mistaken). But again, Oliver is the person to ask.

Best,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

having acquired a somewhat better understanding of how to employ HyperInt in single scale problems, I’m
now trying to learn how to use it when calculating multi-scale integrals.

In principle, many steps are similar to what I did before (using your hints, of course), i.e. usually (assuming linear reducibility)
I can ignore nonlinear terms until the very last integration. When integrating over last the remaining Feynman parameter
the nonlinear terms are clearly important and must be accounted for by adjusting the value of _hyper_splitting_field accordingly.

However, I’ve encountered some cases where I got an intermediate result with one last integration remaining, yet
failed to get to the correct final result from there. Here is one example

expr:=[[1/(x[4]*mm1), [[-(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(-mm2*x[4]-mm2), -(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(mm1*x[4]-mm2*x[4]-mm2)]]], [-1/((x[4]+1)*mm1), [[-x[4]-1, -(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(mm1*x[4]-mm2*x[4]-mm2)]]], [1/(x[4]*mm1), [[-(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(-mm2*x[4]-mm2), -(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(-mm2*x[4]+mm1-mm2)]]], [-1/(x[4]*mm1), [[-(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(-mm2*x[4]-mm2)], [-(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(mm1*(x[4]+1))]]], [1/(x[4]*mm1), [[-(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(mm1*(x[4]+1))], [-x[4]-1]]], [1/(x[4]*mm1), [[0, -(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(-mm2*x[4]-mm2)]]], [-1/((x[4]+1)*x[4]*mm1), [[0, -(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(mm1*x[4]-mm2*x[4]-mm2)]]], [-1/((x[4]+1)*mm1), [[0, -(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(-mm2*x[4]+mm1-mm2)]]], [-1/((x[4]+1)*x[4]*mm1), [[-x[4]-1, -(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(-mm2*x[4]+mm1-mm2)]]], [-1/(x[4]*mm1), [[-(-mm2*x[4]^2+mm1*x[4]-2*mm2*x[4]-mm2)/(-mm2*x[4]-mm2), -x[4]-1]]]]:

When evaluating

forgetAll():
_hyper_splitting_field :={}:
_hyper_abort_on_divergence:=false;
_hyper_ignore_nonlinear_polynomials := false;
_hyper_verbosity:=4;
_hyper_algebraic_roots:=false;
_hyper_verbose_frequency:=20;
forgetAll():
hyperInt(expr,x[4]):

HyperInt complains about a polynomial that doesn’t factor linearly in x[4]. Then I add

_hyper_splitting_field :={RootOf(mm2*_Z^2 + (-mm1 + 2mm2)_Z + mm2)}:

and run the same code again. In this case I get

Error, (in linearFactors) -(mm2*RootOf(mm2*_Z^2+(-mm1+2*mm2)*_Z+mm2)-mm1+2*mm2)/mm2 is not a polynomial, something went wrong in transformWord

If I activate

_hyper_algebraic_roots:=true;

I end up with yet another error

Error, (in reglimWord) cannot determine if this expression is true or false: FAIL < infinity

Naively, I would have expected that if I have only one integration left, it should be possible to express the result in terms of some (possibly ugly) GPLs.
Is there something I’m missing here?

I’m also attaching a more complete example to the very first post. Here mm1 and mm2 are some masses squared. Normally, one is much larger than the other one,
so I could of course expand, yet here I’d like to understand how to tackle the full integral.

Cheers,
Vladyslav

[Vladyslav Shtabovenko]

I guess the last two commits are not related to my issue. In any case activating cgSingleReductionClique does not alter
the error messages.

[Erik Panzer]

Dear Vladyslav,

the 2nd error is thrown because HyperInt does not know how to compute regularized limits of hyperlogarithms with algebraic (non-rational) arguments. In the particular case of your integral, there would be an easy workaround, but I am planning to address this more generally (which will take time).

The 1st error is thrown simply to abort calculations where such roots appear, precisely because that will cause problems in other parts. I’ll probably remove that first error message though, as it is somewhat misleading.

The easiest way around this problem would be if you rationalized the root, which is easy to do in this case.

I understand that it would be useful to remove more error messages such that you can always force a result out of HyperInt for such cases were all but the last integration are rational. I put it on the todo list.

Best,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

sorry for the late reply and many thanks for your hints and code updates. Indeed, I can rationalize my square roots
in mm1 and mm2 as follows:

resAux2Ep0:=eval(factor(eval(resAux1Ep0,[mm1=1/(t[1]*t[2]),mm2=(1-t[2]^2)/(4*t[1]*t[2])]))):

This way I can obtain the O(ep^0) result which nicely agrees with pySecDec (for mm1=1 and allowed values of mm2).
However, when trying to calculate the O(ep^1) piece, I encountered yet another kind of error:

Error, (in reglimWord) Pinch-singularity not implemented; need -delta[t[2]]=delta[t[2]] |HyperInt.mpl:372|

One can reproduce the issue by taking the expression from https://pastebin.com/8Q7vNw26
and evaluating

hyperInt(expr,x[4]):

The full code is also attached (r is sqrt(1-4 mm2/mm1) and I set mm1=1). Is there some workaround you could suggest for this type of errors?

PS Not sure if helpful, but I think that at least for cases with square roots one could employ the RationalizeRoots
package by Besier, Wasser and Weinzierl. In my example the package readily suggests the appropriate transformation

<< RationalizeRoots`
RationalizeRoot[Sqrt[mm1^2 - 4 mm1 mm2]]
% // InputForm

being

{{mm1 -> 1/(t[1]*t[2]), mm2 -> -1/4*(-1 + t[2]^2)/(t[1]*t[2])}}

Cheers,
Vladyslav

[Erik Panzer]

Dear Vladyslav,

I’ll have to investigate this further. In the meantime, I found that changing the parametrization to

[mm1=1,mm2=-t*(1+t)]

seems to work (you can reconstruct the dependence on mm1 by power counting). Does this work for you too?

Best,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

sorry for the late reply. Somehow during all the “Weihnachtstrubel“ I forgot to answer here, although I intended to.

With your transformation I can indeed recover the correct results (cross-checked against pySecDec). It’s interesting that in this
case one cannot apply fibrationBasis with respect to r as in

fibrationBasis(eval(resEp1v1,[t=(r-1)/2]),[r])

since this leads to

Error, (in reglimWord) Pinch-singularity not implemented; need -delta[r]=delta[r] |HyperInt.mpl:372|

On the other hand, one can of course do it the other way around as in

eval(fibrationBasis(resEp1v1,[t]),t=(r-1)/2):

which doesn’t produce any errors.

Although your transformation works, I was also wondering if one can find such transformations in a more systematic way,
since not everyone has a good eye for such variable changes.

Luckily, it seems that the RationalizeRoots package can generate further transformations if one uses the options
MultipleSolutions and ForceFDecomposition. Indeed,

RationalizeRoot[Sqrt[1 - 4 mm2], MultipleSolutions -> True, ForceFDecomposition -> True] 

suggests

{{mm2 -> -1/4*(1 + 2t[1])/t[1]^2}, {mm2 -> -2t[1](1 + 2t[1])}}

with the second transformation being identical to your suggestion (up to a rescaling t[1]->t[1]/2). Needless to say
that with this choice I can also recover the correct result.

To sum it up, it seems that RationalizeRoots can be very useful in conjunction with RationalizeRoots when handling
integrals that depend on several variables. I’m attaching my updated notebooks, in the hope that they might be useful to fellow HyperInt users.

Many thanks for your help. It didn’t occur to me that one might need to choose between different variable transformations
to find the one that works best with HyperInt.

‌

[Erik Panzer]

Dear Vladyslav,

thanks for confirming that this works! And also thanks for uploading your RationlizeRoots files, I hope they will be helpful to others.

Many thanks for your help. It didn’t occur to me that one might need to choose between different variable transformations
to find the one that works best with HyperInt.

In principle this should not matter and any linearizing variable change should do (and indeed is good enough for integration). In fibrationBasis, however, the program has to pick consistent branches, and for some parametrizations the origin (r=0) can correspond to pinch-singularities. You should always get rid of these problems by a simple translation of the origin (which is what t vs r does). Hopefully, I’ll get at some point to implement also the pinch case in full generality, to avoid this annoyance.

‌

With best wishes,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

let me report another glitch that I encountered while using HyperInt.

Here I'm performing the very last integration of a 3-loop 7-line Feynman parameter
integral to obtain the O(ep)-piece. Quite surprisingly, this fails
after about 60 seconds with the following error message

Error, (in fibrationBasis) assigning to a long list, please use Arrays |HyperInt.mpl:1710|

To my understanding this looks like some programming issue, rather than
a poor behavior of the integral itself.

I'm attaching a minimal working example to reproduce the problem, but
if needed I could also provide a full worksheet that leads to auxResEp(1).

Cheers,
Vladyslav

[Vladyslav Shtabovenko]

  <div class="preview-container wiki-content"><!-- loaded via ajax --></div>
  <div class="mask"></div>
</div>

</div> </form>

[Vladyslav Shtabovenko]

OK, the bug seems to be related to the mechanism of storing/checking for the divergences. Setting

_hyper_check_divergences:=false:

before doing the integration I can get the final result just fine. The divergence check is not so relevant
in my case, as I anyhow cross-check all the results numerically.

[Vladyslav Shtabovenko]

And here is another issue I found during my HyperInt sessions. In the attached example,
when using the order of integrations suggested by suggestIntegrationOrder, I run into

Error, (in convert/parfrac) denominator factors must contain the parfrac variable |HyperInt.mpl:680|

during the O(ep) integration in the last variable. However, if I switch two last variables in my
vars list (contrary to what suggestIntegrationOrder tells me), the integration works just fine
and delivers the correct result.

‌

[Vladyslav Shtabovenko]

attached example

[Erik Panzer]

Thanks for reporting! I fixed the “Error, (in fibrationBasis) assigning to a long list, please use Arrays |HyperInt.mpl:1710|" in <<cset 6fde05b98b92>>

[Erik Panzer]

Your other problem (partfracBug.mw) is due to the ill-definedness of the integral. Already at earlier stages there are divergences, and you get i*pi*delta[x[k]]'s. I could change the code to ignore such terms, but that would be dangerous, as it could suggest to users that everything is fine, while actually the integral is ill-defined, so I prefer the current abort that is triggered by the error.

If you are in a situation where you can make sense of these ill-defined quantities, and you want to force the integration through, you can do the following:

fibrationBasis(aux,[x[4]]):
eval(%,delta[x[4]]=del4):
res:=hyperInt(%,x[4]);

The renaming of delta[x[4]] to del4 masks it from hyperInt, so it treats del4 just as some parameter, instead of getting worried of the indeterminacy during the integration of x[4]. Does this solve your issue?

Overall, my feeling is that the proper way out of this is to define the integrals and imaginary parts such that no divergences occur and everything is well-defined without cuts inside the integration domain (in particular the positive real axis, which is the integration domain). In your example, aux jumps between the lower and the upper half-plane, so it’s not clear which branch an integral along the real axis should pick.

‌

[Vladyslav Shtabovenko]

Many thanks for the explanation. In principle, I’m all in favor of doing things the proper way,
but it’s not always clear how to do so when using a particular code.

So, what would be the standard procedure to define the imaginary parts
in the Feynman parameter representation such, that HyperInt would know which branch to take?
Perhaps there are some examples I could look at?

[Erik Panzer]

In the documentation it is explained what the delta[x] mean (they pick the half-plane in case of cuts on the positive real axis). But it is not clear if you can pick the contour you want by specifying those branches. In particular, it may be that Feynman i*eps prescription does not correspond to a contour with fixed delta’s.

In your concrete example, the situation looks bad - there are a lot of minus signs in your 2nd Symanzik polynomial (F), and depending on which variables are small or large, the imaginary part of F will be dictated by different Schwinger/Feynman parameters. Hence the leading monomial, determining the imaginary part, switches throughout the integration domain, and sometimes it has positive, sometimes negative sign. I could not immediately see a pattern among those two classes of monomials that would fix the imaginary part (via delta’s) in a way that it has constant sign everywhere.

Probably the easiest way to avoid this problem is by introducing more parameters (masses, momenta) until you have a non-empty Euclidean region.

[Vladyslav Shtabovenko]

I understand the source of delta[x], but my question is rather if one can pick a contour “by hand“ within HyperInt.
If I would be calculating something similar in Mathematica, which by default always selects +I eta for logs or square
roots, I would simply add an “+/-I*eta” (with eta being an undefined symbol) to the FP representation. This effectively prevents
Mathematica from taking the branch cut in quantities such as log(-x) so that after expanding in ep one can select the branch
cut by hand. This is how I usually handle this problem for integrals that are simple enough to be expressed as a product of
Gamma functions.

In the case of HyperInt, I guess that adding an extra symbol times an “I” would mess up the internal integration routines, so it’s probably
not the right way to proceed. On the other hand, tools using sector decomposition usually can get the sign of the imaginary
part right without extra tricks, so perhaps I should look into the corresponding codes in more details.

Introducing more parameters surely would help, but I’m afraid that it will also make the integration much more complicated
if not impossible (e.g. if an integral is “simple” when on-shell but becomes elliptic off the mass shell). Or perhaps I misunderstood
your suggestion?

My current approach of “fitting” the delta’s from numerical results is admittedly not a well-defined way to calculate Feynman
integrals, but in practice it seems to work quite well (at least for the integrals I’m currently dealing with). At least at the moment
I’m not aware of another procedure that would allow me to get analytic results with as little effort as I have now (this is why
I’m now a huge fan of HyperInt and advertise it to all my colleagues as a very powerful tool to tackle linearly reducible Feynman integrals).
But I’m always eager to learn how to do things properly.

Concerning the error, I’m just wondering whether the error message
Error, (in convert/parfrac) denominator factors must contain the parfrac variable |HyperInt.mpl:680|
is clear enough to an average user. As you explained, the underlying cause is an issue with the imaginary parts, but
to me the error rather seems to suggest that it’s something about a nonlinearity in the denominators or so. Perhaps it would more
descriptive to write that (in this case) it’s because of the delta’s appearing inside functions that need to be integrated.

I fully agree with you that HyperInt shouldn’t ignore such cases, hence pretending that everything is fine, while it’s not, but on the other
hand it’s also obvious that one cannot always detect such issues. In my example, switching the two last integration variables makes
the integration go through, although the final result is still in-defined in your language.
May be one could consider adding an option that would ignore the appearance of the delta’s in such situations,
hence catering for cases where the user knows how to make sense of potentially ill-defined results? Of course one would
add all due warnings in the description of such an option. But this is certainly only up to you to decide. The workaround you
posted is sufficient for my purposes and I’m also very grateful for the explanation of the underlying issue.

[Vladyslav Shtabovenko]

Dear Erik,

in my attempts to explore what is possible with HyperInt and what is not, I’ve got a question that hopefully
should be of interest also for other HyperInt users.

From my current (admittedly limited) understanding, if one gets the “Error, (in partialFractions) XYZ is not linear in XXX |HyperInt.mpl:673|”
error message despite of having _hyper_ignore_nonlinear_polynomials set to true, this signals that some denominators in the arguments of Hlogs
contain polynomials nonlinear in the current integration variable, so that HyperInt cannot carry out the integration. This may actually occur
only during intermediate integrations, since in the very last integration (when there is only integration variable left), one can always employ
_hyper_splitting_field:=RootOf(...):to get through.

Now, having such nonlinearities that raise the above mentioned error may signal that the integral in question is not linearly reducible. In that
case one clearly cannot use HyperInt to calculate it. However, as far as I understand, in some cases the integral might still be linearly reducible
and one just needs to find a proper variable transformation to linearize the “bad“ denominators. There is a talk by R. Schabinger where he did
this in a very nontrivial case, but the integrals I’m having in mind are, in fact, much simpler than that.

To have a specific example, consider e.g. a 2-loop on-shell integral with 3 massive lines forming a triangle. This single scale integral has been calculated many years ago,
and can be found e.g. in https://arxiv.org/pdf/hep-ph/9907431.pdf where the integral is denoted as F10101 (cf. Fig. 1). It is finite
and the O(ep^0) result is given on p.18 of the linked paper. The results is also rather simple, containing only zeta(2), zeta(3) and a Clausen
function. So the O(ep^0) piece should be certainly calculable via direct Feynman parametrization (or am I wrong here?)

Trying to recover this result with HyperInt I can first derive the FP-representation of this integral

(x[4]*(x[3] + x[5]) + x[1]*(x[2] + x[3] + x[5]) + x[2]*(x[3] + x[4] + x[5]))^(-1 + 3*ep)*(x[1]^2*(x[2] + x[3] + x[5]) + x[2]^2*(x[3] + x[4] + x[5]) + x[4]*(x[3]*(x[4] - x[5]) + x[4]*x[5]) + 
   x[1]*(x[2]^2 + x[3]*(x[4] - x[5]) + x[4]*x[5] + x[2]*(2*x[3] + x[4] + x[5])) + x[2]*(x[3]*(x[4] - x[5]) + x[4]*(x[4] + 2*x[5])))^(-1 - 2*ep)

Since it is not projective, I projectivize it in the simplest way via x[i]->x[i]/(x[1]+…x[5]) . Trying to naively calculate this integral with HyperInt one actually gets stuck on the
second to last integration (see attachment), where the infamous “Error, (in partialFractions) XYZ is not linear in XXX |HyperInt.mpl:673|” pops up.

So my question is, how one should approach such cases when using HyperInt? The integral is clearly a very simple one, but HyperInt seemingly gets stuck here.
Should one then try to find a suitable integration variable transformation that linearizes the offending denominators, or are there other way one could circumvent the problem?

[Erik Panzer]

Dear Vladyslav,

it is usually indeed the case that, if you get

“Error, (in partialFractions) XYZ is not linear in XXX |HyperInt.mpl:673|”
error message despite of having _hyper_ignore_nonlinear_polynomials set to true, this signals that some denominators contain polynomials nonlinear in the current integration variable, so that HyperInt cannot carry out the integration.

Note that this is about the denominators of the rational functions multiplying hyperlogarithms, not about arguments of the hyperlogarithms. It can happen that non-linear denominators appear which are spurious (cancel between various terms in a big expression). But in most cases, as you say, indeed this is a sign that the Landau variety of the integrand is not linear in the integration variable.

Concretely, in your example:

The short answer is that your integral is not linearly reducible in your parametrization, so if you want to integrate it with HyperInt, you have to change variables.

By the way, I do not understand what you mean by

Since it is not projective, I projectivize it in the simplest way via x[i]->x[i]/(x[1]+…x[5]) . 

The integrand you gave is already projective!

But going on: I set ep=0 so the integrand is 1/(U*F) where (I just rewrite x[i] as xi to get rid of all those square brackets)

U:=x1*x2+x1*x3+x1*x5+x2*x3+x2*x4+x2*x5+x3*x4+x4*x5;
F:=(x2+x3+x5)*(x1^2+x4^2)+(x1+x3+x4+x5)*x2^2
+(x1*x2-x1*x3+x1*x4-x2*x3+2*x2*x4-x3*x4)*x5+2*x1*x2*x3+x1*x2*x4+x1*x3*x4+x2*x3*
x4;

The polynomial reduction finds 2 linear integrations:

L:=table(): S:=irreducibles({U,F}):
L[{}]:=[S,combinat[choose](S,2)]:
cgReduction(L):
suggestIntegrationOrder(L);

I get a suggestion of [x3,x5]. In the reduction

L[{x3,x5}][1];

you see 4 polynomials that have quadratic terms:

Q:={x1^2+x1*x2+x1*x4+x2*x4+x4^2, 
x1^2+x1*x2+x1*x4+x2^2+2*x2*x4+x4^2, 
x1^2+2*x1*x2+x1*x4+x2^2+x2*x4+x4^2, 
x1^2+2*x1*x2+x1*x4+x2^2+2*x2*x4+x4^2};

This is not linearly reducible, because these polynomials do not all factor by merely extending the base field. For example, say I want to integrate x2 next, then we get square roots of the discriminants

map(discrim, Q, x2);

which gives

{1, -3*x1^2, -3*x4^2, 4*x1*x4}

So by adjoining _hyper_splitting_field:={I,sqrt(3)}, you can take all the roots, except the last one: sqrt(4*x1*x4)=2*sqrt(x1*x4) remains algebraic.

To rationalize, let’s set x1=1 and substitute x4=x4^2, then this root becomes 2*sqrt(x4^2)=2*x4. So all roots are rational, hence you can integrate x2 after this substitution, and in the end you can integrate x4 (for that integration, I had to set _hyper_abort_on_divergence:=false, because HyperInt could not verify symbolically that there is a cancellation - but you can verify by looking into _hyper_divergencesthat the expression is indeed zero, so the integral is convergent).

As a result this gives

1.56883858 - I*Pi^3/18*delta[x4]

where the real part is a long expression (I do not have a reduction table for 3rd roots of unity at hand). Numerically, this seems to agree perfectly with the expression for F10101 on page 18 in the paper that you linked, provided you integrate x4 with a negative imaginary part (below the real axis) so that delta[x4]=-1.

‌

Does this help?

Best,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

many thanks for your detailed answer!

First of all, I agree that the integral is already projective: I had a bug in my code for checking the projectivity property,
which I discovered few days after my posting. But then I didn't want to create additional confusion/presure by follow up comments,
so I left the statement as it is. In the meantime, I think that the most convenient way to check projectivity is to substitute the
Feynman parameters variables by some large random primes times the rescaling parameter. This seems to be much less error-prone
than trying to show that the rescaling parameter cancels out algebraically by using Simplify and alike. This "RandomPrime-approach" is currently implemented here.

If someone is interested, my workflow is to use FeynCalc with

int = FAD[p2]*FAD[p2 + q]*FAD[{p1, m1}]FAD[{p1 - p2, m1}]FAD[{p1 + q, m1}]

aux = FCFeynmanParametrize[int, {p1, p2}, Names -> x, Indexed -> True,FCReplaceD -> {D -> 4 - 2 ep}, Simplify -> True,Assumptions -> {ep > 0},FinalSubstitutions -> {FCI@SPD[q] -> mm, m1^2 -> mm,m2^2 -> mmc}] /. mm | mmc -> 1

{psi, phi} =FCFeynmanPrepare[int, {p1, p2}, Names -> x, Indexed -> True,FinalSubstitutions -> {FCI@SPD[q] -> mm, m1^2 -> mm,m2^2 -> mmc}][[1 ;; 2]] /. mm | mmc -> 1

auxP = FCFeynmanProjectivize[aux[[1]], x, {}]

After that I evaluate

psi // InputForm
phi // InputForm
auxP // InputForm

and paste the output into my Maple session with HyperInt

Regarding the generation of L[{}], the code L[{}]:=[S,combinatchoose]:indeed seems to do a much better job than just L[{}] := [S, {S}];that is given in most of the examples that come with HyperInt. I realized that when trying to recalculate
some 2-loop on-shell integrals (J001 and J111) from the aforementioned paper and failing already during the cgReduction(L)-part. Going through the HyperInt documentation pdf-file I found your hint on p.26, which then did the job for me. For F10101 there seems to be no difference, but I suppose that it is safe to assume that L[{}]:=[S,combinatchoose]is a better default choice?

Regarding my original question concerning the treatment of such integrals that do not work out of the box, recently there has been
a new preprint dealing precisely with this issue and explaining some neat tricks:

https://arxiv.org/pdf/2103.15423.pdf

But as your kind explanation shows, most "simpler" integrals essentially require only a change of variables that rationalizes the roots
of the corresponding polynomials. I still have a feeling that this procedure could be automatized (especially using the RationalizeRoots package), but perhaps I need to work out more examples first, until I can see a clear pattern here.

For the sake of other users I collected the code necessary to calculate F10101 according to your guide in the attached notebook.

One strange thing I noticed is that the integration in x[4] seems to fail at first with

Error, (in reglimWord) cannot determine if this expression is true or false: min(0, FAIL) < FAIL |HyperInt.mpl:333|

However, if I simply repeat the previous integration in x[2] and then evaluate the code to integrate in x[4] again, then miraculously it works. Perhaps you could have a look at it and share your thoughts on this behavior.

Another thing I wanted to ask, is whether one could add an automatic conversion of HyperInt zeta's to the Maple format when using evalf? Something similar to my zetaToZeta helper function so that one doesn't have to do it by hand. That would make the numerical checks of the HyperInt output much more convenient.

[Erik Panzer]

Dear Vladyslav,

that’s a neat idea to check projectivity. Indeed unfortunately it can be tricky to get this reliably detected symbolically.

Automatization of changing variables is a tricky business. What I had been pondering for some time is adding a routine to the polynomial reduction that is triggered when non-linear polynomials are encountered. In the case of only quadrics, it would compute all discriminants, and could then report those, which would already be helpful (in fact, this is what I always do ‘by hand’, following the steps I explained in my previous post).

I would however refrain from implementing automatic changes of variables, because if this is done repeatedly, and ‘bad’ choices are made, the expression size can explode. Also, in some cases, one can find a scaling like xi->xi*xj or such for some variables that renders a set of polynomials linear in a further variable. Such a transformation would be much preferred over a quadratic reparametrization of discriminants. Overall, it is an unsolved and ill-defined problem to find ‘the best’ parametrization.

Another thing I wanted to ask, is whether one could add an automatic conversion of HyperInt zeta's to the Maple format when using evalf?

I think the best solution would be to replace the indices zeta[…] by parantheses zeta(…). This way, it would be possible to extend evalf of Maple to recognize zeta(…) and thereby completely dispose of the need for any conversion. (Conversion to Zeta(…) does not work in general, since Zeta(n,m) in Maple is already used, but not for the double zeta value - rather the mth derivative of the Riemann (single) zeta).

Indeed, if I would rewrite HyperInt now, I would certainly go for zeta(…) instead of zeta[…] I think. At the time, I did not know about the subtleties of the Maple language that make zeta[…] less flexible here.

Unfortunately this is a problem with backward compatibility. If I change zeta[…] to zeta(…), then old code relying on zeta[…] notation won’t run with new HyperInt’s.

I totally agree that it would be nice to have zeta’s seamlessly integrated with evalf, but due to the notation problem I do not yet know how to achieve that without breaking something.

Best,

Erik

[Erik Panzer]

Hello again,

regarding

Error, (in reglimWord) cannot determine if this expression is true or false: min(0, FAIL) < FAIL |HyperInt.mpl:333|

This is caused because you used _hyper_algebraic_roots, so that you got Root() expressions in the integration of x[2], as you see in the output of resAux2Ep0. HyperInt cannot work with any expressions containing Root()s. When you recomputed the x[2] integral, in your file you had also set

_hyper_splitting_field:={I,sqrt(3)}:

The effect is that instead of Root()'s, the expression of the resAux2Ep0 now is in terms of I and sqrt(3). This is an expression that HyperInt understands, so the subsequent integration over x[4] can proceed.

As a rule of thumb I suggest:

Never use _hyper_algebraic_roots:=true; except in the last integration.

I agree that the current error message is pretty meaningless. I should improve it and make HyperInt complain explicitly about Root’s in the input.

Best,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

many thanks for the clarifications. I agree that it may not be tractable to automatize the variable change in full generality,
but perhaps one could aim at some special cases?

Anyhow, it would be nice to have more explicit examples dealing with particular cases, so that one can try to apply the
methods shown there to the problem at hand. Perhaps something like a HyperInt example gallery
(a pun to https://feyncalc.github.io/examples). If you are interested, I would be happy to contribute some examples
reproducing literature results. In any case, I have a strong feeling that HyperInt is very much underestimated in the HEP
community: there are few people who employ it brilliantly to tackle tough problems (I don’t belong to the chosen ones),
some people who have heard of it but never tried it (or they erroneously think that it’s not useful for what they do) and a lot
of people who have no idea that it exists. And this is a shame, IMHO.

Regarding the zeta’s, you are of course right that it’s somewhat tricky. For example, even though Maple now supports MZVs,
all the indices must be positive. So to evaluate/simplify something like zeta[1,-3] I currently need to export it to PolyLogTools in
Mathematica (where it is denoted as MZV[{-3,1}]). Are there some plans to include such constants via additional tables or so?
I guess that one could scrap them from the HPL package (which is what PolyLogTools relies upon) and convert the output to
the Maple format.

Cheers,
Vladyslav

‌

[Erik Panzer]

Dear Vladyslav,

Perhaps something like a HyperInt example gallery   (a pun to https://feyncalc.github.io/examples). If you are interested, I would be happy to contribute some examples reproducing literature results.

That’s a great idea indeed! I agree it could be very useful to have such a gallery. The reason that it does not exist already is partly that I do not know how to technically manage and set this up in a practical and sustainable way:

• Worksheet examples are the most immediately useful within Maple. However, they may cause problems for people with older versions of Maple. Also, they are not nice to view outside Maple or online.
• Plain text Maple code is my go-to way to use HyperInt, but that may put off some people. It needs good documentation, or ideally needs to be explained in snippets step-by-step.
• It would be nice to have the examples online, not just as a download of a file. Would it make sense to have them as wiki entries, or standalone webpages, or rendered markdowns?
• Another approach would be to add an appendix to the documentation, where each example gets a section on its own. That way, all would be in one place, viewable online, and one could easily mix code with explanatory text.

If you have any thoughts/experience/suggestions on these issues, that would be very helpful! I think the format/workflow needs to be fixed in advance, before it makes sense to put time into cleaning up code and making examples presentable and useful.

So to evaluate/simplify something like zeta[1,-3] I currently need to export it to PolyLogTools in Mathematica (where it is denoted as MZV[{-3,1}]).

The zeta[…] notation does not work with evalf due to the parantheses vs bracket issue that I mentioned. But you can write them as multiple polylogarithms explicitly:

evalf(Mpl([1,3],[1,-1]));
# output: .8778567139e-1-.1986033885e-9*I

Does this work for you?

With best wishes,

Erik

[Vladyslav Shtabovenko]

Dear Erik,

Worksheet examples are the most immediately useful within Maple. However, they may cause problems for people with older versions of Maple. Also, they are not nice to view outside Maple or online.

Plain text Maple code is my go-to way to use HyperInt, but that may put off some people. It needs good documentation, or ideally needs to be explained in snippets step-by-step.

I agree. In Mathematica there is the same issue with plain text .m-files vs fancy notebooks in .nb-format and I think that one should keep only plain text scripts in the repository.
I think that in the case of HyperInt the idea would be to have all the examples as plain text. For generating the gallery they could be loaded into Maple and exported as HTML files.

It would be nice to have the examples online, not just as a download of a file. Would it make sense to have them as wiki entries, or standalone webpages, or rendered markdowns?

I would say that Markdown is the most convenient thing because it nicely fits into static websites that one can generate on BitBucket

https://support.atlassian.com/bitbucket-cloud/docs/publishing-a-website-on-bitbucket-cloud/

Maple actually offers the possibility to export evaluated worksheets to LaTeX. From there one could use pandoc to generate Markdown files and add them to a static website.

Another approach would be to add an appendix to the documentation, where each example gets a section on its own. That way, all would be in one place, viewable online, and one could easily mix code with explanatory text.

In principle yes, but PDF files are not so well indexed by search engines as websites. One could also do both once there is a set of suitable examples available.

If you have any thoughts/experience/suggestions on these issues, that would be very helpful! I think the format/workflow needs to be fixed in advance, before it makes sense to put time into cleaning up code and making examples presentable and useful.

One should certain think of some kind of template that will be used to treat the example integrals. 90% of steps are always the same, so that having a template can significantly facilitate the process of
adding new examples.

‌

The zeta[…] notation does not work with evalf due to the parantheses vs bracket issue that I mentioned. But you can write them as multiple polylogarithms explicitly:

evalf(Mpl([1,3],[1,-1]));
# output: .8778567139e-1-.1986033885e-9*I

Many thanks! Indeed, I didn’t know that this works. But in this case I don’t really understand why one can’t have an automatic evalf working with the output of HyperInt.
Is this only about the brackets? Couldn’t one do something like

zetaToAuxZeta:=proc(x):
subsindets(x,'specindex(anything,zeta)',f1->auxZeta(op(f1)))
end proc:
####
zetaToAuxZeta(c1*zeta[1,2]+c2*zeta[2])

to get rid of the square brackets and then convert all auxZeta functions to appropriate functions that Maple can evaluate numerically (e.g. Zeta and your Mpl). What am I missing here?

Cheers,
Vladyslav

‌

‌

‌

‌

‌