\documentclass[12pt]{article}
\usepackage[nosort]{cite}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{amsmath}
%\usepackage{axodraw}
%\usepackage{epsfig}
%\usepackage{varioref}
\usepackage{latexsym}
\usepackage{amssymb}
\usepackage{amsfonts}
\textwidth 6.0in
\textheight 8.5in
\topmargin -0.25truein
\oddsidemargin 0.30truein
\raggedbottom
\def\e{{\rm e}}
\newcommand{\E}{\mbox{e}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\al}{\alpha}
\newcommand{\Be}{\mbox{B}}
\newcommand{\gm}{\gamma}
\newcommand{\Gm}{\Gamma}
\newcommand{\dl}{\delta}
\newcommand{\Dl}{\Delta}
\newcommand{\ep}{\varepsilon}
\newcommand{\sg}{\sigma}
\newcommand{\kp}{\kappa}
\newcommand{\ka}{\kappa}
\newcommand{\lm}{\lambda}
\newcommand{\Lm}{\Lambda}
\newcommand{\om}{\omega}
\newcommand{\pa}{\partial}
\newcommand{\dd}{\mbox{d}}
\newcommand{\dr}{{\rm d}}
\newcommand{\la}{\langle}
\newcommand{\ra}{\rangle}
\newcommand{\MS}{\mbox{MS}}
\newcommand{\nn}{\nonumber}
\newcommand{\D}{\mbox{d}}
\newcommand{\lra}{\leftrightarrow}
\newcommand{\bfm}[1]{\mbox{\boldmath$#1$}}
\newcommand{\I}{i}
\newcommand{\insl}{\not\!\, \in}
\newcommand{\nisl}{\not\!\, \ni}
\newcommand{\cdo}{\!\cdot\!}
\newcommand{\ut}{\underline{t}}
\newcommand{\uta}{\underline{\tau}}
\newcommand{\ual}{\underline{\alpha}}

\begin{document}
\parindent=1.5pc

\begin{flushright}
  \begin{tabular}{@{}l@{}}
  TTK-12-21 \\
  TTP12-015 \\
  SFB/CPP-12-31
  \end{tabular}
\end{flushright}

\vspace{0.5cm}

\begin{center}
  \textbf{\LARGE Expansion by regions:
    revealing potential \\[1ex] and Glauber regions automatically}

  \vspace{1cm}

  \textbf{\renewcommand{\thefootnote}{\alph{footnote}}%\large
    Bernd Jantzen$^{1,}$\footnote{e-mail: jantzen@physik.rwth-aachen.de},
    Alexander V. Smirnov$^{2,4,}$\footnote{e-mail: asmirnov80@gmail.com} and
    Vladimir A. Smirnov$^{3,4,}$\footnote{e-mail: smirnov@theory.sinp.msu.ru}
    \setcounter{footnote}{0}
  }

\vspace{-0.5cm}

\textit{%\small
\begin{list}{}{\setlength{\leftmargin}{0.5cm}\setlength{\itemindent}{-\leftmargin}
  \setlength{\topsep}{0cm}\setlength{\itemsep}{0cm}}
\item $^1$Institut f\"ur Theoretische Teilchenphysik und Kosmologie,
  RWTH Aachen University, 52056 Aachen, Germany
\item $^2$Scientific Research Computing Center,
  Moscow State University, 119992 Moscow, Russia
\item $^3$Skobeltsyn Institute of Nuclear Physics, Moscow State University,
  119992 Moscow, Russia
\item $^4$Institut f\"ur Theoretische Teilchenphysik, KIT,
  76128 Karlsruhe, Germany
\end{list}}

\vspace{0.5cm}

\begin{abstract}
  \noindent
  When performing asymptotic expansions using the strategy of expansion by
  regions, it is a non-trivial task to find the relevant regions. The
  recently published Mathematica code \texttt{asy.m} automates this task,
  but it has not been able to detect potential regions in threshold
  expansions or Glauber regions. In this work we present an algorithm and
  its implementation in the update \texttt{asy2.m} which also reveals
  potential and Glauber regions automatically.
\end{abstract}

\end{center}

%\vspace{0.5cm}



\section{Introduction}

If a given Feynman integral depends on kinematic invariants and masses
which essentially differ in scale, a natural idea is to expand it
in ratios of small and large parameters. As a result, the integral is written
as a series of simpler quantities than the original integral itself and it can be
substituted by a sufficiently large number of terms of such an expansion.
For limits typical of Euclidean space (for example, the off-shell large-momentum limit
or the large-mass limit), one can write down the corresponding asymptotic expansion
in terms of a sum over certain subgraphs of a given
graph~\cite{Chetyrkin:1988zz,Chetyrkin:1988cu,Gorishnii:1989dd,
Smirnov:1990rz,Smirnov:1994tg,books1a}. This prescription of expansion by subgraphs
has been mathematically proven (see \cite{Smirnov:1990rz} and Appendix~B.2 of
\cite{books1a}).
Moreover, there is an automated tool~\cite{Seidensticker:1999bb,Harlander:1997zb}
where such an expansion by subgraphs is implemented.

For limits typical of Minkowski space (i.e.\ which cannot be formulated in
Euclidean space) the universal strategy of
expansion by regions~\cite{Beneke:1997zp,Smirnov:1998vk,Smirnov:1999bza,books1a}
is available. It consists of the following prescriptions:
\begin{itemize}
\item
  Divide the space of the loop momenta into various regions and, in
  every region, expand the integrand in a Taylor series with respect
  to the parameters that are considered small there.
\item
  Integrate the integrand, expanded in the appropriate way in every
  region, over the \emph{whole integration domain} of the loop momenta.
\item
  Set to zero any scaleless integral.
\end{itemize}
%
As shown in~\cite{Smirnov:1999bza}, this prescription can also be applied
to parametric representations of Feynman integrals, i.e.\ alpha parameters
(or generalized Feynman parameters) integrated from 0 to~$\infty$,
eventually restricted by a delta function. Then the regions are specified
by scaling relations between the parameters.

There is no mathematical proof that this prescription is correct in all
situations. But also no examples are known where a proper application of
the expansion by regions leads to wrong results. An indirect proof exists
for limits typical of Euclidean space because here the strategy of regions
is equivalent to the mathematically proven expansion by subgraphs.
%
A systematic study of the expansion by regions was presented recently
in~\cite{Jantzen:2011nz}. There it was shown explicitly and
illustrated using various one-loop examples that one can start from a
decomposition of a given integral into non-intersecting domains and arrive
at an expansion by regions in the above sense. This requires certain
conditions on the choice and completeness of the considered regions which
were derived in~\cite{Jantzen:2011nz}. As pointed out there, the appearance
of additional overlap contributions can be avoided by adequate choices of
the regions and regularization parameters.

While these findings provide some hints on the proper choice of the
regions, it remains a non-trivial task to actually reveal the typical
regions for a given limit. Usually, one starts from considering one-loop
examples, checks the results against known analytical results, then
proceeds in two loops etc.  One can also use the second version
\cite{FIESTA2} of the code {\tt FIESTA} \cite{FIESTA} to obtain numerically
several first terms of a given asymptotic expansion.

Recently an algorithm for an automatic search of regions was suggested and
implemented on a computer as the open source {\tt Mathematica} code
\texttt{asy.m}~\cite{asy}.
%
The algorithm uses a geometric approach based on finding the convex hull
of a set of points determined from a parametric representation of the
Feynman integral. In this way all possible sets of scalings for the
(Feynman) parameters are found which lead to non-vanishing (because
non-scaleless) integrals. These regions may then be used to expand the
parametric integral, or they can be translated into regions for expanding
the integral in loop-momentum space.

This code works successfully for a large variety of limits, at least in
cases where the function~$\mathcal{F}$ in the corresponding parametric
representation which depends on kinematic invariants and masses is only
composed of terms with the same sign. Moreover, it was shown in \cite{asy}
that in this case there are no regions except for the ones produced by the
code (all others result in scaleless integrals). In particular, the code
{\tt asy.m} works for Sudakov-type limits which are typical of Minkowski
space. As it was pointed out in \cite{asy} the code does not reveal
potential regions in threshold expansions, and similarly it fails to detect
the so-called Glauber regions.

Our first goal in the present paper is to provide an algorithm and the
corresponding update {\tt asy2.m} which automatically reveals potential
regions. This is done in Section~\ref{sec:potential}.
%
In Section~\ref{sec:Glauber} we consider the problem of revealing Glauber
regions. Here the situation is more complicated than for the potential
regions, but the updated code {\tt asy2.m} is now also able to reveal
Glauber regions automatically.
%
In both Sections \ref{sec:potential} and~\ref{sec:Glauber}, we first
formulate simple changes of variables and decompositions of a given Feynman
integral, using instructive one-loop examples, which lead to integrals
where {\tt asy2.m} is able to detect the relevant regions corresponding to
scalings of the parameters. Then, for both cases, we explain how to use the
new features of {\tt asy2.m} to perform these algorithmic steps
automatically.

In the case of Section~\ref{sec:Glauber} with Glauber regions, the
structure of the regions differs depending on whether the expansion is
performed in loop-momentum space or at the level of the parametric
integral. We show in Section~\ref{sec:Glauber_l} how to disentangle and
match the various regions arising in this problem using generic propagator
powers, and how {\tt asy2.m} can be used to automate such an analysis.

A summary of the new features and the syntax of \texttt{asy2.m} (together
with a download link) is provided in Section~\ref{sec:syntax}. In
Section~\ref{sec:conclusion} we conclude by discussing the mathematical
problem of proving the expansion by regions for a simple example which is
not related to Feynman integrals, but where \texttt{asy2.m} works
successfully.



\section{Revealing potential contributions}
\label{sec:potential}

Let us consider the one-loop propagator diagram with two massive lines in
the threshold limit, i.e.\ when $y=m^2-q^2/4\to 0$ with $q$ being the external
momentum:
\bea
F(q^2,m^2)=\int \frac{\dd^d k}{(k^2-m^2)\; \bigl((k-q)^2-m^2\bigr)}\;,
\label{prop1_integral}
\eea
where $d=4-2\ep$, and the usual $+i 0$ is implied in all the propagators.
Within the strategy of expansion by regions, the hard and the potential
regions give contributions to the expansion \cite{Beneke:1997zp,books1a}.
The previous version of the code {\tt asy.m} \cite{asy} reported only about
the hard region. The reason for this can be seen in the corresponding
parametric representation,
\bea
F(q^2,y)=i\pi^{d/2} \, \Gm(\ep)
%&&\nn \\  && \hspace*{-50mm} \times
\iint \frac{(x_1+x_2)^{2\ep-2}\;\delta\left(x_1+x_2-1\right)\;
\dd x_1\dd x_2
}{\left[
\frac{q^2}{4}(x_1-x_2)^2 + y(x_1+x_2)^2
-i 0\right]^{\ep}}\;,
\label{prop1_alpha}
\eea
where the parameters~$x_i$ are integrated from 0 to $\infty$ (restricted by
the delta function).
As it was pointed out in \cite{asy}, it is the region where
$x_1\approx x_2$ (more precisely \mbox{$x_1-x_2 \sim y^{1/2}$}) which causes problems.
In other words, the polynomial in the square brackets
in~(\ref{prop1_alpha}) (considered at positive $q^2$ and $y$) has terms of
different sign, such that cancellations occur because of the presence
of the negative term $-q^2 x_1x_2/2$.

To reveal the missing potential contribution, let us perform a simple trick.
We decompose the integration domain into two subdomains, $x_1\leq x_2$ and
$x_2\leq x_1$. The two resulting integrals are equal to each other, but
such an equality will not generally take place for any integral.
In the first domain we turn to new variables by
$x_1=x_1'/2,\;x_2=x'_2+x'_1/2$, remove the primes at $x_i$ and obtain
the integral (again from 0 to $\infty$ with the usual restrictions
via the delta function)
\bea
i\pi^{d/2} \, \frac{\Gm(\ep)}{2}
\iint \frac{(x_1+x_2)^{2\ep-2}\;\delta\left(x_1+x_2-1\right)\;
\dd x_1\dd x_2}
{\left[
\frac{q^2}{4}x_2^2 + y(x_1+x_2)^2
-i 0\right]^{\ep}}\;.
\label{prop1_alpha_1}
\eea
The goal of this trick was to make the line $x_1=x_2$ (in the old variables)
the border of an integration domain which turned out to be
(in the new variables) $x_2=0$. Now we can run the code {\tt asy2.m}.
Since this is a parametrical integral rather than a Feynman integral
we use the newly introduced command\footnote{The name of the command refers
  to its application to parametric integrals contributing to Wilson loops.}
{\tt WilsonExpand[]} for integrals where all parameters are integrated from
0 to $\infty$:
\begin{verbatim}
  WilsonExpand[q^2/4*x2^2 + y*(x1 + x2)^2, x1 + x2,
    {x1, x2}, {q -> 1, y -> x}, Delta -> True]
\end{verbatim}
The first two arguments of \texttt{WilsonExpand[]} are the polynomials
$\mathcal{F}$ and $\mathcal{U}$ from the parametric
integral~(\ref{prop1_alpha_1}) as defined in~\cite{asy}. The third argument
is the list of integration parameters, and the fourth argument specifies
the scaling of the kinematic quantities with respect to the small parameter
which is labelled by the global symbol~\texttt{x}. The option
\verb@Delta -> True@ tells \texttt{WilsonExpand[]} that, under the
integration, the sum over an arbitrary non-empty subset of the integration
parameters is restricted to~1 by a delta function.

Note that \texttt{WilsonExpand[]} can only take into account such a delta
function if the specific choice of the sum over parameters in the argument
of the delta function is irrelevant. This is the case for the generalized
Feynman-parameter integrals (\ref{prop1_alpha}) and~(\ref{prop1_alpha_1}),
where we could e.g.\ replace $\delta(x_1+x_2-1)$ by $\delta(x_1-1)$ without
changing the value of the integral. But if a specific form of the delta
function is assumed, e.g.\ by replacing $x_1+x_2 \to 1$ under the integral,
then the option \texttt{Delta} of \texttt{WilsonExpand[]} does not apply
(see Section~\ref{sec:syntax} for details). Alternatively,
\texttt{WilsonExpand[]} can be used without the option \texttt{Delta} after
eliminating one of the integrations, e.g.\ via $\delta(x_1-1)$.

The output of \texttt{WilsonExpand[]} is a list of regions specified by the
scaling of the parameters~$x_i$. Here the result \verb@{{0, 0}, {0, 1/2}}@
reveals not only the hard region with the scaling $x_1,x_2\sim y^0$ but
also the region $x_1\sim y^0,x_2\sim y^{1/2}$ which we qualify as the
potential region.

The contribution of the potential region can then be evaluated
in terms of gamma functions at general $\ep$. Taking into account that
we have two identical integrals after our decomposition, we arrive at
the following result for the potential contribution:
\be
i \pi^{d/2} \, \Gm(\ep-1/2)
\sqrt{\frac{\pi y}{q^2}} \, y^{-\ep}
\label{prop1_pot}
\ee
in agreement with~\cite{books1a}.

In fact, such a trick of making manifest squares of some linear combination
of the integration parameters was already used in the code
{\tt FIESTA} \cite{FIESTA} in order to evaluate numerically Feynman integrals
at a threshold. Using the implementation of this procedure in {\tt FIESTA}
it turned out to be possible to automate the above trick for a general Feynman
integral.
%
In the present version {\tt asy2.m}, the user may call the command
\texttt{AlphaRepExpand[]} with the additional option \texttt{PreResolve}
enabled which automatically looks for the change of variables described above:
\begin{verbatim}
  AlphaRepExpand[{k}, {k^2 - m^2, (k - q)^2 - m^2},
    {q^2 -> qq, m^2 -> qq/4 + y}, {qq -> 1, y -> x},
    PreResolve -> True]
\end{verbatim}
As in the previous version of \texttt{asy.m}, the arguments of
\texttt{AlphaRepExpand[]} are the list of loop momenta, the list of
denominators of the loop integral, a list of replacements for the kinematic
quantities, and the list of scalings with respect to the small
parameter~\texttt{x}.
The output is a list of entries which indicate for each region the changes of
variables, the Jacobian of the integral transformation, and the scalings of
the new variables:
\begin{verbatim}
  {{{x[1] -> y[1] + y[2]/2, x[2] -> y[2]/2}, 2, {0, -1/2}},
   {{x[1] -> y[1] + y[2]/2, x[2] -> y[2]/2}, 2, {0, 0}},
   {{x[1] -> y[1]/2, x[2] -> y[1]/2 + y[2]}, 2, {0, 0}},
   {{x[1] -> y[1]/2, x[2] -> y[1]/2 + y[2]}, 2, {0, 1/2}}}
\end{verbatim}
Remember that a region remains invariant if all parameters are rescaled
simultaneously, due to the homogeneity of the polynomials $\mathcal{F}$
and~$\mathcal{U}$. So the scaling $\{0,-1/2\}$ shown above is equivalent to
$\{1/2,0\}$, i.e.\ the first new parameter is suppressed by $y^{1/2}$ with
respect to the second one.
%
For both parts of the integral decomposition the hard and potential regions
are found.

The preresolution algorithm implemented in {\tt asy2.m} (and switched on
with the option \texttt{PreResolve}) tries to eliminate factorized
combinations of terms in the function $\mathcal{F}$ which potentially
cancel each other, like $(x_1-x_2)^2$ in the example above. It checks all
pairs of variables (say, $x$ and $y$) which are part of monomials with
opposite sign. For all those pairs the code tries to build a linear
combination $z$ of $x$ and $y$ such that in the variables $x$ and $z$ or
$y$ and $z$ this monomial disappears. The code checks whether in the new
variables the number of monomials with opposite sign decreases. For all
such pairs the code recursively repeats the initial procedure in the new
variables. As a result it creates a tree of possible bisections and
corresponding replacements of variables. A leaf of this tree is a set of
sectors and functions such that one cannot decrease the number of monomials
with opposite sign any longer. Ideally it means that all monomials now have
the same sign. The code analyzes all leafs and chooses one of those with
the minimal number of opposite-sign monomials (or the minimal number of
sectors if the numbers of monomials with opposite sign coincide).
%
After finishing with the preresolution, the code performs the replacements
and looks for regions in all those sectors.

Note that the algorithm can only find the necessary variable
transformations if it is able to determine the relative signs of all terms
in the polynomial~$\mathcal{F}$. As the signs of symbols are unknown to
\texttt{asy2.m}, the substitution rules of the third and fourth arguments
of \texttt{AlphaRepExpand[]} must replace all kinematic quantities by
numbers (integers or fractions of integers) or powers of the small
parameter~\texttt{x}.

We have checked that the updated version {\tt asy2.m} works in various
examples of the threshold expansion (considered in~\cite{Beneke:1997zp,books1a}):
a triangle, a box, the two-loop
propagator diagram (with the masses $m,m,m,m,0$), a two-loop vertex
diagram. Because of the decomposition of a given integration domain into
subdomains, the number of resulting integrals for various regions increases a little
bit. For example, the (hard-hard) region for the two-loop
propagator diagram is described by six integrals, the (potential-ultrasoft)
region is also described by six integrals, etc. However, the
(potential-hard) region is described by four integrals with some regions
(with scalings composed of 1, 1/2 and 0), and four more integrals with a set
of regions of a different type (composed of 1 and 0).

Let us finally mention that the preresolution algorithm also works for
threshold expansions with unequal masses. Returning to the one-loop
example~(\ref{prop1_integral}), but now with different masses $m_1$ and
$m_2$ in the propagators and the expansion parameter defined as $y =
(m_1+m_2)^2/4 - q^2/4 \to 0$, we call:
\begin{verbatim}
  AlphaRepExpand[{k}, {k^2 - m1^2, (k - q)^2 - m2^2},
    {q^2 -> (m1 + m2)^2 - 4*y},
    {y -> x, m1 -> 2, m2 -> 3/2},
    PreResolve -> True]
\end{verbatim}
We have set the values of the masses to different rational numbers, $m_1
\to 2$ and $m_2 \to 3/2$ (the actual values are irrelevant), thus
permitting the preresolution algorithm to know about their sign without
assuming any equality or other relation between them. The
polynomial~$\mathcal{F}$ in the parametric representation of the Feynman
integral reads $\mathcal{F} = (m_1 x_1 - m_2 x_2)^2 + 4y x_1 x_2$, and
cancellations occur where $m_1 x_1 = m_2 x_2$. These cancellations are
automatically made explicit by adequate changes of variables. The output
\begin{verbatim}
  {{{x[1] -> y[1] + 3*y[2]/7, x[2] -> 4*y[2]/7}, 7/4, {0, -1/2}},
   {{x[1] -> y[1] + 3*y[2]/7, x[2] -> 4*y[2]/7}, 7/4, {0, 0}},
   {{x[1] -> 3*y[1]/7, x[2] -> 4*y[1]/7 + y[2]}, 7/3, {0, 0}},
   {{x[1] -> 3*y[1]/7, x[2] -> 4*y[1]/7 + y[2]}, 7/3, {0, 1/2}}}
\end{verbatim}
shows that \texttt{asy2.m} detects both regions. The variable
transformations are always normalized such that sums of the parameters
(here $x_1+x_2$) remain invariant.

So, we now have a manifestly Lorentz-invariant treatment of threshold
expansion and a code that automatically provides the set of relevant regions.



\section{Revealing Glauber contributions}
\label{sec:Glauber}

\begin{figure}[t]
  \centering
  \includegraphics[scale=0.9]{Glauber-5point}
  \caption{One-loop five-point integral exhibiting a Glauber contribution.}
  \label{fig:Glauber}
\end{figure}
Let us consider the one-loop five-point integral in Fig.~\ref{fig:Glauber},
where two initial-state partons both perform a collinear splitting into two
partons each with momenta $p_1,p_2$ and $q_1,q_2$, respectively. While two
partons, one of each pair, collide with a large centre-of-mass energy
$Q=\sqrt{(p_2+q_2)^2}$, the two remaining partons exchange a particle with
the small mass~$m$.
We will use the simplified kinematics
$p_1 = p_2 = p$ and $q_1 = q_2 = q$ with $p^2 = q^2 = 0$
and $(p+q)^2 = 2p\cdot q = Q^2$ in the limit  $m^2/Q^2 \to 0$:
\bea
F(Q^2,m^2)=\int \frac{\dd^d k}{(k^2-m^2)(k^2-2p\cdot k) (k^2+2p\cdot k)
(k^2-2q\cdot k) (k^2+2q \cdot k)}\;.
\label{5pt_integral}
\eea

Before we search for regions using \texttt{asy.m}, we notice that this
five-point integral is similar to the Sudakov form factor example
treated in Section~6 of~\cite{Jantzen:2011nz}. From the viewpoint of the
convergence of the expansions, the second and third propagators
of~(\ref{5pt_integral}) are equivalent, and so are the fourth and fifth
propagators. Effectively, the five-point integral has only three different
types of propagators which are equivalent to the ones of the three-point
integral in~\cite{Jantzen:2011nz}. So the integral~(\ref{5pt_integral}) can
be expanded in loop-momentum space employing the regions and convergence
domains known from~\cite{Jantzen:2011nz} (and using generic propagator
powers as analytic regulators where necessary):
\begin{itemize}
\item a hard region where $k \sim Q$,
\item a 1-collinear region where $k^2 \sim p\cdot k \sim m^2$ and $q\cdot k
  \sim Q^2$,
\item a 2-collinear region where $k^2 \sim q\cdot k \sim m^2$ and $p\cdot k
  \sim Q^2$,
\item a Glauber region where $p\cdot k \sim q\cdot k \sim m^2$, and the
  components of~$k$ perpendicular to the plane spanned by $p,q$ scale as
  $k_\perp \sim m$.
\end{itemize}
The collinear-plane region mentioned in~\cite{Jantzen:2011nz} yields
only scaleless contributions. But, in contrast to the three-point integral,
the five-point integral has a non-vanishing Glauber contribution. The
Glauber region even provides the leading contribution scaling as
$(m^2)^{-2-\ep}$, whereas the collinear contributions start with
$(m^2)^{-1-\ep}$ and the hard contribution starts with $(m^2)^0$.

The five-point integral~(\ref{5pt_integral}) can be represented in terms of
an integral over Feynman parameters,
\begin{align}
F(Q^2,m^2)&=-i\pi^{d/2} \, \Gm(3+\ep)
\nn \\ &\times
\idotsint \frac{(x_1+\ldots+x_5)^{1+2\ep}\;\delta\left(\sum_i x_i-1\right)\;
\dd x_1\cdots\dd x_5
}{\bigl[x_1(x_1+\ldots+x_5) m^2
+(x_2-x_3)(x_4-x_5)Q^2-i 0\bigr]^{3+\ep}
 }\;,
\label{Fp_integral}
\end{align}
where one can choose the sum in the argument of the delta function in an
appropriate way, i.e.\ restrict only the sum over a subset of the
parameters to 1 and extend the integration over the rest of the parameters
to the whole domain $[0,\infty)$.

Applying the strategy of expansion by regions in Feynman-parameter space
and trying to reveal regions relevant to the given limit $m^2/Q^2 \to 0$
with the help of the code {\tt asy.m} \cite{asy}, we call:
\begin{verbatim}
  AlphaRepExpand[{k},
    {k^2 - m^2, k^2 - 2*p*k, k^2 + 2*p*k, k^2 - 2*q*k, k^2 + 2*q*k},
    {p^2 -> 0, q^2 -> 0, p*q -> Q^2/2},
    {Q -> 1, m^2 -> x}]
\end{verbatim}
The output states the following set of three regions:
\begin{verbatim}
  {{0, 0, 0, 0, 0}, {0, 0, 0, 1, 1}, {0, 1, 1, 0, 0}}
\end{verbatim}
As before, the regions are specified by the scaling of the Feynman parameters
in terms of powers of the small parameter~$m^2$. For example, for the
second region we have
$x_1 \sim x_2 \sim x_3 \sim(m^2)^0$, $x_4 \sim x_5 \sim(m^2)^1$.
The first region is hard; its contribution starts with $(m^2)^0$.
The second and third regions start with order $(m^2)^{-1-\ep}$. They
correspond to the two collinear regions stated for the momentum-space
expansion above. But \texttt{asy.m} does not find anything corresponding to
the Glauber region; in particular, none of the regions found by
\texttt{asy.m} provides the leading $(m^2)^{-2-\ep}$ contribution.

We notice that, as in the previous section about potential contributions,
the polynomial in the square brackets of~(\ref{Fp_integral}) has terms of
different sign. The missing Glauber contribution stems from the parameter
region where either $(x_2-x_3) \sim (m^2)^1$ or $(x_4-x_5) \sim (m^2)^1$. So let us
decompose the parametric integral into four parts corresponding to the
domains where the two factors $(x_2-x_3)$ and $(x_4-x_5)$ are either
positive or negative and then introduce new variables in such a way that
this product takes the form $\pm x'_2 x'_4$. For example, in the domain
$x_2\leq x_3$, $x_5\leq x_4$ we change the variables by
$x_2=x'_3/2,\;x_3=x'_2+x'_3/2$ and by $x_4=x'_4+x'_5/2,\;x_5=x'_5/2$,
similarly to our example in the previous section.
%
However, in the threshold expansion the cancelling terms appeared in
squared form such that a transformation between one pair of variables was
sufficient. Here two separate factors involve cancellations, which requires
a twofold change of variables.

Removing the primes from the variables~$x_i$, the parametric integral reads
\bea
  F(Q^2,m^2)=2(I_+ + I_-)
\eea
with
\begin{align}
I_{\pm} &=-i\pi^{d/2} \, \frac{\Gm(3+\ep)}{4}
\nn \\  &\times
\idotsint \frac{(x_1+x_2+x_3+x_4+x_5)^{1+2\ep}\;\delta\left(x_1-1\right)\;
\dd x_1\cdots\dd x_5
}{\left[x_1(x_1+x_2+x_3+x_4+x_5) m^2
\pm x_2 x_4Q^2 -i 0\right]^{3+\ep}
 }\;,
\label{Fp_integral_A}
\end{align}
where we have chosen the argument of the delta function as $x_1-1$, so that
we may also write
\begin{align}
I_{\pm} &=-i\pi^{d/2} \, \frac{\Gm(3+\ep)}{4}
\nn \\  &\times
\int_0^\infty \cdots \int_0^\infty
\frac{(1+x_2+x_3+x_4+x_5)^{1+2\ep}\;
\dd x_2\cdots\dd x_5
}{\left[(1+x_2+x_3+x_4+x_5) m^2
\pm x_2 x_4Q^2 -i 0\right]^{3+\ep}
 }\;.
\label{Fp_integral_B}
\end{align}
It is sufficient to consider the expansion of $I_+$ and obtain a result for
$I_-$ by analytically continuing $Q^2 \to -Q^2-i0$, taking into account
that the dependence on $Q^2$ is power-like.

Now we can apply {\tt asy2.m} to the integral~$I_+$ using either the command
\begin{verbatim}
  WilsonExpand[x1*(x1 + x2 + x3 + x4 + x5)*m^2 + x2*x4*Q^2,
    x1 + x2 + x3 + x4 + x5,
    {x1, x2, x3, x4, x5},
    {Q^2 -> 1, m^2 -> x}, Delta -> True]
\end{verbatim}
for integrals~(\ref{Fp_integral_A}) restricted by a delta function (see
Section~\ref{sec:syntax}), yielding the output:
\begin{verbatim}
  {{0, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 0, 0, 1, 0}}
\end{verbatim}
Or we use the command
\begin{verbatim}
  WilsonExpand[(1 + x2 + x3 + x4 + x5)*m^2 + x2*x4*Q^2,
    1 + x2 + x3 + x4 + x5,
    {x2, x3, x4, x5},
    {Q^2 -> 1, m^2 -> x}]
\end{verbatim}
for integrals~(\ref{Fp_integral_B}) over variables from $0$ to $\infty$
without any restriction, and obtain the output:
\begin{verbatim}
  {{0, 0, 0, 0}, {1, 0, 0, 0}, {0, 0, 1, 0}}
\end{verbatim}
The two results are equivalent, as $x_1 \sim (m^2)^0$ is implied in the
second case. So we obtain again three regions. We will see in a moment that
this list of regions is indeed correct and complete.

But first, let us emphasize that the new preresolution algorithm in
\texttt{asy2.m} is capable of performing the transformation of the integral
from (\ref{Fp_integral}) to~(\ref{Fp_integral_A}) automatically:
\begin{verbatim}
  AlphaRepExpand[{k},
    {k^2 - m^2, k^2 - 2*p*k, k^2 + 2*p*k, k^2 - 2*q*k, k^2 + 2*q*k},
    {p^2 -> 0, q^2 -> 0, p*q -> Q^2/2},
    {Q -> 1, m^2 -> x},
    PreResolve -> True]
\end{verbatim}
The output of this command lists four different variable transformations
according to the twofold decomposition described above. For each of the
integrals over new parameters, the regions $\{0,0,0,0,0\}$, $\{0,1,0,0,0\}$
and $\{0,0,0,1,0\}$ are found, up to permutations in the order of the
parameters from different changes of variables.

The evaluation of the contributions to each region, as found by
\texttt{WilsonExpand[]} or \texttt{AlphaRepExpand[]} (including the
\texttt{PreResolve} option), is straightforward. The first region is the
hard one. The contributions of the second and third regions are not
individually regularized by dimensional regularization, as it often happens
for Sudakov-type limits. We use an auxiliary analytic regularization by
introducing additional powers $x_2^{\dl_2} x_3^{\dl_3} x_4^{\dl_4}
x_5^{\dl_5}$ of the new variables into the integrand
of~(\ref{Fp_integral_B}), taking the limit $\dl_2,\dl_3,\dl_4,\dl_5 \to 0$
in the end. The leading-order (LO) contribution of the second and third
regions to the integral $F(Q^2,m^2)$ reads
\begin{align}
  -i\pi^{d/2} \, \frac{i\pi \Gamma(\ep)}{2Q^2(m^2)^{2+\ep}}\;.
  \label{Fp_LO}
\end{align}
This agrees with the leading contribution of the Glauber region in the
momentum-space expansion.

We have found the leading Glauber contribution of order
$(m^2)^{-2-\ep}$. But we seem to have lost the two collinear regions with
the scalings $\{0,0,0,1,1\}$ and $\{0,1,1,0,0\}$ found before the change of
variables. In fact, we can evaluate the contributions from these two
regions by expanding the integral~(\ref{Fp_integral_B}). The resulting
integrals are scaleless and regularized by the parameters $\dl_3, \dl_5$,
so they vanish, and \texttt{asy2.m} is right in omitting these two regions.

We are also able to solve the integral~(\ref{Fp_integral_B}) including the
auxiliary analytic regularization factor $x_2^{\dl_2} x_3^{\dl_3}
x_4^{\dl_4} x_5^{\dl_5}$ in terms of a onefold Mellin--Barnes
representation:
\begin{align}
I_{\pm} &=-i\pi^{d/2} \, \frac{\Gm(1+\dl_3) \Gm(1+\dl_5)}{4} \,
\frac{1}{2\pi i} \int \dd z \, (m^2)^z (\pm Q^2-i0)^{-3-\ep-z}
\nn \\  &\times
\Gm(-z) \Gm(-2-\ep+\dl_2-z) \Gm(-2-\ep+\dl_4-z)
\nn \\  &\times
\frac{\Gm(1-\dl_2-\dl_3-\dl_4-\dl_5+z) \Gm(3+\ep+z)}{\Gm(-1-2\ep-z)}
\;.
\label{Fp_integral_MB}
\end{align}
The asymptotic expansion of~$I_\pm$ in the limit $m^2/Q^2 \to 0$ is
obtained by taking the residues of the poles of the functions
$\Gm(\ldots-z)$. The poles of $\Gm(-z)$ correspond to the hard region,
while the poles of the two functions $\Gm(-2-\ep+\dl_{2,4}-z)$ provide the
contributions of the second and third regions. So \texttt{asy2.m} has found
all contributing regions.

In the Mellin--Barnes integral~(\ref{Fp_integral_MB}) we can safely take
the limit $\dl_2,\dl_3,\dl_4,\dl_5 \to 0$, add up $I_+$ and $I_-$, and
arrive at the Mellin--Barnes representation
\begin{align}
F(Q^2,m^2) &=i\pi^{d/2} \, \frac{i}{2} \, \frac{1}{2\pi i}
\int \dd z \, (m^2)^z (Q^2)^{-3-\ep-z} e^{i\pi(\ep+z)/2}
\nn \\  &\times
\Gm(-z) \Gm(-2-\ep-z) \Gm\left(\tfrac{-1-\ep-z}{2}\right) \,
\frac{\Gm(1+z) \Gm\left(\frac{3+\ep+z}{2}\right)}{\Gm(-1-2\ep-z)}
\;.
\label{Fp_integral_MB2}
\end{align}
The LO contribution to $F(Q^2,m^2)$ is obtained from the residue of the
single pole at $z=-2-\ep$, in agreement with~(\ref{Fp_LO}).
%
The next-to-leading-order (NLO) contribution stems from the residue of the double
pole at $z=-1-\ep$ and reads
\bea
i\pi^{d/2} \, \frac{ \Gamma (1+\ep)}{(Q^2)^2(m^2)^{1+\ep}}
\left(i \frac{\pi}{2}
+2 \psi(-\ep)- \psi(1+\ep)+ \gm_{\rm E}-\ln\frac{Q^2}{m^2} -1\right)\;.
\label{Fp_NLO}
\eea
This agrees with the NLO contributions of the second and third regions.

At next-to-next-to-leading order (NNLO) there is a contribution from the
residue of the single pole at $z=-\ep$ which reads
\bea
-i\pi^{d/2} \, \frac{i\pi \Gamma (2+\ep)}{4(Q^2)^3(m^2)^{\ep}}
\label{Fp_NNLO_g}
\eea
and agrees with the NNLO contributions of the second and third regions. The
second NNLO contribution comes from the residue of the single pole at
$z=0$. It is given by
\bea
  -i\pi^{d/2} \,
  \frac{i \, e^{i\pi\ep/2} \, \Gm(-2-\ep) \Gm\left(\frac{1+\ep}{2}\right)
      \Gm\left(\frac{1-\ep}{2}\right)}{
    2 (Q^2)^{3+\ep} \, \Gm(-1-2\ep)}
\label{Fp_NNLO_h}
\eea
and agrees with the LO contribution of the hard region.
%
So indeed all contributions to the five-point integral up to NNLO are
correctly reproduced by the contributions of the three regions found by
\texttt{asy2.m} after the decomposition of the integral and the change of
variables.

When revealing Glauber regions for a general diagram, the preresolution
algorithm of \texttt{asy2.m} tries to eliminate monomials with opposite
sign in the polynomial~$\mathcal{F}$ by automatically separating the
integration into domains and performing changes of variables. If the option
\texttt{PreResolve} is enabled for \texttt{AlphaRepExpand[]}, the code
warns the user if the elimination of monomials with opposite sign has not
been successful, such that possibly not all regions are revealed. This is
the case if some monomials of opposite sign remain in the
polynomial~$\mathcal{F}$ after tries to eliminate them, or if symbols with
unknown signs are present in the polynomial. We are therefore convinced
that \texttt{AlphaRepExpand[]}, with the preresolution enabled, either
reveals all relevant regions or issues a warning.

As for the threshold expansion of the previous section, \texttt{asy2.m} can
also treat more complicated kinematical situations, e.g. the five-point
integral depicted in Fig.~\ref{fig:Glauber} with $p_1 \ne p_2$ and $q_1 \ne
q_2$ (retaining $p_1 \parallel p_2$ and $q_1 \parallel q_2$ such that $p_i
\cdot p_j = q_i \cdot q_j = 0$ and $2p_i \cdot q_j = P_i Q_j$ with
$P_i,Q_j>0$):
\begin{verbatim}
  AlphaRepExpand[{k},
    {k^2 - m^2, k^2 - 2*p1*k, k^2 + 2*p2*k,
                k^2 - 2*q2*k, k^2 + 2*q1*k},
    {p1^2 -> 0, p2^2 -> 0, p1*p2 -> 0,
      q1^2 -> 0, q2^2 -> 0, q1*q2 -> 0,
      p1*q1 -> P1*Q1/2, p1*q2 -> P1*Q2/2,
      p2*q1 -> P2*Q1/2, p2*q2 -> P2*Q2/2},
    {m^2 -> x, P1 -> 1, P2 -> 3, Q1 -> 2, Q2 -> 3/2},
    PreResolve -> True]
\end{verbatim}
All kinematic invariants are replaced by rational numbers in order to
enable the preresolution algorithm to work. The code correctly decomposes
the integral into four pieces and finds the three regions for each of them.

Let us finally discuss the reason why, besides the hard region which is always
present, the expansion in loop-momentum space requires two collinear
regions and one Glauber region, whereas the expansion of the parametric
integrals~(\ref{Fp_integral_B}) has two regions providing the leading
Glauber contribution and no further collinear regions.
%
In fact, the momentum-space expansion is also valid for loop
integrals~(\ref{5pt_integral}) where each propagator is raised to an
arbitrary, even non-integer power. For the decomposition of the parametric
integral and the change of variables, however, we have assumed the specific
form~(\ref{5pt_integral}) with each propagator present exactly once. In
this case, the loop integrand can be expanded into partial fractions as follows:
\begin{align}
  \lefteqn{\frac{1}{(k^2-m^2)(k^2-2p\cdot k) (k^2+2p\cdot k)
    (k^2-2q\cdot k) (k^2+2q \cdot k)}}
  \qquad
\nn \\ &
  = \frac{1}{4 (m^2)^2}
    \left( \frac{1}{k^2-m^2} - \frac{1}{k^2} - \frac{m^2}{(k^2)^2} \right)
    \left( \frac{1}{k^2-2p\cdot k} + \frac{1}{k^2+2p\cdot k} \right)
\nn \\ &\times
    \left( \frac{1}{k^2-2q\cdot k} + \frac{1}{k^2+2q\cdot k} \right)
  \;.
\end{align}
Expanding this product of terms, one obtains twelve three-point integrals,
which are well known. Because they only depend on $Q^2 = 2p\cdot q$, we
recognize from the last two factors the structure $F(Q^2,m^2) = 2(I_+ +
I_-)$, where $I_+$ and $I_-$ are related by $Q^2 \to -Q^2$ as before. The
three-point integrals with the massless propagators $1/k^2$ or $1/(k^2)^2$
only have a hard region. The massive three-point integral with propagator
$1/(k^2-m^2)$ is known to possess a hard and two collinear regions. Its LO
and NLO hard contributions are cancelled by the massless three-point
integrals, such that the uncancelled hard contributions start with
$(m^2)^2/(m^2)^2 = (m^2)^0$, as for the five-point integral. The LO
collinear contributions of the three-point integrals, enhanced by the
$1/(m^2)^2$ prefactor, scale as $(m^2)^{-2-\ep}$.

So, in the special case when all propagator powers are equal to~1, the
five-point integral reduces to a linear combination of three-point
integrals revealing the same structure of regions as found from the
expansion of the parametric integral~(\ref{Fp_integral_B}).
This picture changes when generic propagator powers are introduced as
analytic regulators, which is done in the next section.



\section{Disentangling regions via propagator powers}
\label{sec:Glauber_l}

In the previous section we have seen different patterns of regions
arising when expanding either in loop-momentum space or in parametric
space. The individual contributions can be disentangled more easily when
the dependence on the propagator powers is retained. Instead
of~(\ref{5pt_integral}), let us consider the integral
\begin{align}
  F(Q^2,m^2)&=\int \frac{\dd^d k}{(k^2-m^2)^{1+\lm_1}
    (k^2-2p\cdot k)^{1+\lm_2} (k^2+2p\cdot k)^{1+\lm_3}}
\nn \\ &\times \frac{1}{
    (k^2-2q\cdot k)^{1+\lm_4} (k^2+2q\cdot k)^{1+\lm_5}}\;,
  \label{5pt_integral_l}
\end{align}
where the analytic regularization parameters~$\lm_i$ make the propagator
powers different from the previous case. The asymptotic expansion in
loop-momentum space yields contributions from the four regions listed in
the beginning of Section~\ref{sec:Glauber}. The LO hard contribution still
scales as $(m^2)^0$, but now the LO 1-collinear contribution scales as
$(m^2)^{-1-\ep-\lm_1-\lm_2-\lm_3}$, the LO 2-collinear contribution as
$(m^2)^{-1-\ep-\lm_1-\lm_4-\lm_5}$ and the LO Glauber contribution as
$(m^2)^{-2-\ep-\lm_1-\lm_2-\lm_3-\lm_4-\lm_5}$. So all regions are
characterized by a distinct scaling. When performing the expansion in
parametric space for generic~$\lm_i$, we are able to disentangle the
contributions from each region and match them with the regions obtained in
loop-momentum space.

For generic~$\lm_1,\ldots,\lm_5$, the parametric integral corresponding
to~(\ref{Fp_integral}) reads
\begin{align}
\lefteqn{
F(Q^2,m^2)=-i\pi^{d/2} \, e^{-i\pi(\lm_1+\ldots+\lm_5)} \,
  \frac{\Gm(3+\ep+\lm_1+\ldots+\lm_5)}{\Gm(1+\lm_1) \cdots \Gm(1+\lm_5)}
} \quad
\nn \\ &\times
\idotsint \frac{x_1^{\lm_1} \cdots x_5^{\lm_5} \;
(x_1+\ldots+x_5)^{1+2\ep+\lm_1+\ldots+\lm_5}\;
\delta\left(\sum_i x_i-1\right)\;\dd x_1\cdots\dd x_5
}{\bigl[x_1(x_1+\ldots+x_5) m^2
+(x_2-x_3)(x_4-x_5)Q^2-i 0\bigr]^{3+\ep+\lm_1+\ldots+\lm_5}
 }\;.
\label{Fp_integral_l}
\end{align}
The change of variables performed in Section~\ref{sec:Glauber}, e.g.\ $x_2
= x'_2 + x'_3/2$, $x_3 = x'_3/2$ for $x_2 \ge x_3$, is complicated by the
presence of the factors $x_2^{\lm_2} x_3^{\lm_3} x_4^{\lm_4} x_5^{\lm_5}$,
where parts of the monomials will change into polynomials. If we want to
keep the simple form $F(Q^2,m^2) = 2(I_+ + I_-)$ from the previous section,
then we have to require $\lm_3 = \lm_2$ and $\lm_5=\lm_4$. Under this
restriction of the parameters~$\lm_i$, the parametric integrals can be
written as
\begin{align}
I_\pm&=-i\pi^{d/2} \, \frac{e^{-i\pi(\lm_1+2\lm_2+2\lm_4)}}{4^{1+\lm_2+\lm_4}} \,
  \frac{\Gm(3+\ep+\lm_1+2\lm_2+2\lm_4)}{
    \Gm(1+\lm_1) \Gm^2(1+\lm_2) \Gm^2(1+\lm_4)}
\nn \\ &\times
\idotsint \frac{x_1^{\lm_1} (2x_2+x_3)^{\lm_2} x_3^{\lm_2}
  (2x_4+x_5)^{\lm_4} x_5^{\lm_4} \;
(x_1+\ldots+x_5)^{1+2\ep+\lm_1+2\lm_2+2\lm_4}
}{\bigl[x_1(x_1+\ldots+x_5) m^2
\pm x_2 x_4 Q^2-i 0\bigr]^{3+\ep+\lm_1+2\lm_2+2\lm_4}
 }
\nn \\ &\qquad\qquad{}\times
\delta\left(x_1-1\right)\;\dd x_1\cdots\dd x_5
\;.
\label{Fp_integral_l_A}
\end{align}
In order to find the regions for the asymptotic expansion
of~(\ref{Fp_integral_l_A}), we have to provide the additional polynomial
factors $(2x_2+x_3)$ and $(2x_4+x_5)$ to \texttt{WilsonExpand[]}. We can do
this by multiplying the new polynomials to the second argument of the
command:
\begin{verbatim}
  WilsonExpand[x1*(x1 + x2 + x3 + x4 + x5)*m^2 + x2*x4*Q^2,
    (x1 + x2 + x3 + x4 + x5)*(2*x2 + x3)*(2*x4 + x5),
    {x1, x2, x3, x4, x5},
    {Q^2 -> 1, m^2 -> x}, Delta -> True]
\end{verbatim}
The output of this command reads:
\begin{verbatim}
  {{0, 0, 0, 1, 0}, {0, 0, 0, 1, 1},
   {0, 1, 0, 0, 0}, {0, 1, 1, 0, 0},
   {0, 0, 0, 0, 0}}
\end{verbatim}
In addition to the regions present in the analysis of
Section~\ref{sec:Glauber}, we retrieve the two collinear regions with
scalings $\{0,0,0,1,1\}$ and $\{0,1,1,0,0\}$.
%
The updated code \texttt{asy2.m} is capable of taking into account generic
propagator powers automatically when told through the additional option
\texttt{GenericPowers}:
\begin{verbatim}
  AlphaRepExpand[{k},
    {k^2 - m^2, k^2 - 2*p*k, k^2 + 2*p*k, k^2 - 2*q*k, k^2 + 2*q*k},
    {p^2 -> 0, q^2 -> 0, p*q -> Q^2/2},
    {Q -> 1, m^2 -> x},
    PreResolve -> True, GenericPowers -> True]
\end{verbatim}
When the option \texttt{GenericPowers} is enabled, the
polynomial~$\mathcal{U}$ obtained from the loop integral is multiplied by
the product of all Feynman parameters, $x_1 x_2 \cdots$, before the changes
of variables are performed. Some of these additional factors then turn into
polynomials through replacements in the preresolution algorithm, while
others remain monomials and are therefore irrelevant for the analysis of
\texttt{asy2.m}. The call of \texttt{AlphaRepExpand[]} stated above
correctly yields all five regions for each of the variable transformations.

For the evaluation of the Glauber contributions with scalings
$\{0,0,0,1,0\}$ and $\{0,1,0,0,0\}$, an additional analytic regularization
is needed, and we choose to multiply the integrand
of~(\ref{Fp_integral_l_A}) by $x_2^{\dl_2} x_4^{\dl_4}$. (The parameters
$\dl_3,\dl_5$ from Section~\ref{sec:Glauber} are not needed here due to the
presence of $\lm_2,\lm_4$.) The two Glauber contributions are individually
singular in the limit $\dl_2,\dl_4 \to 0$, but this singularity cancels in
the sum of the two contributions.

The leading contribution to the integral $F(Q^2,m^2)$ originates from the
sum of the LO Glauber contributions. It reads
\begin{multline}
  -i\pi^{d/2} \,
  \frac{i \, e^{-i\pi(\lm_1+2\lm_2+2\lm_4)} \,
      \Gm\left(\tfrac{1}{2}+\lm_2\right) \Gm\left(\tfrac{1}{2}+\lm_4\right)
      \Gm(2+\ep+\lm_1+2\lm_2+2\lm_4)}{
    2 Q^2 \, (m^2)^{2+\ep+\lm_1+2\lm_2+2\lm_4} \,
      \Gm(1+\lm_1) \Gm(1+\lm_2) \Gm(1+\lm_4) \Gm(1-\ep)}
\\ \times
  \Gm(-1-\ep-2\lm_2-2\lm_4)
  \;.
  \label{Fp_l_LO}
\end{multline}
This agrees with the LO contribution from the one Glauber region in the
momentum-space expansion of~(\ref{5pt_integral_l}) for $\lm_3 = \lm_2$ and
$\lm_5=\lm_4$. Here, in the expansion of the parametric integrals, we have
two regions producing the Glauber contribution. This is possible because
the contributions from both regions have the same scaling (for
$\dl_2=\dl_4=0$), starting with $(m^2)^{-2-\ep-\lm_1-2\lm_2-2\lm_4}$.
%
For $\lm_i=0$, the result~(\ref{Fp_l_LO}) reproduces~(\ref{Fp_LO}).

Among the NLO contributions to $F(Q^2,m^2)$, we expect NLO Glauber
contributions and LO collinear contributions. However, the NLO Glauber
contributions vanish exactly for general $\lm_1$, $\lm_2=\lm_3$ and
$\lm_4=\lm_5$, due to non-trivial cancellations between the pieces which
contribute to the NLO expansion of~(\ref{Fp_integral_l_A}) for either of
the scalings $\{0,0,0,1,0\}$ and $\{0,1,0,0,0\}$. The same happens in
loop-momentum space, where the NLO Glauber contribution is proportional to
$(\lm_3-\lm_2)(\lm_5-\lm_4)$, thus vanishing in the case considered here.
%
So the NLO contribution to $F(Q^2,m^2)$ is made up entirely from the LO
collinear contributions. The 1-collinear region provides
\begin{multline}
  -i\pi^{d/2} \,
  \frac{e^{-i\pi(\lm_1+2\lm_2+2\lm_4)} \, e^{i\pi\lm_4} \,
      \Gm(\lm_2-\lm_4) \Gm(1-2\lm_4)
      \Gm(1+\ep+\lm_1+2\lm_2)}{
    2 (Q^2)^{2+2\lm_4} \, (m^2)^{1+\ep+\lm_1+2\lm_2} \,
      \Gm(1+\lm_1) \Gm(1+\lm_2) \Gm(1-\lm_4)}
\\ \times
  \frac{\Gm(-\ep-2\lm_2)}{\Gm(-\ep-2\lm_4)} \;
  \frac{1}{1+2\lm_4}
  \;,
  \label{Fp_l_NLO_1c}
\end{multline}
in agreement with the momentum-space expansion. The 2-collinear
contribution is obtained from this by exchanging $\lm_2 \leftrightarrow
\lm_4$. Adding the two collinear contributions together and performing the
limit $\lm_1,\lm_2,\lm_4 \to 0$, the result~(\ref{Fp_NLO}) is reproduced.

Considering finally the NNLO contributions to $F(Q^2,m^2)$, we expect NNLO
Glauber contributions, NLO collinear contributions and a LO hard
contribution. But here the NLO collinear contributions vanish exactly for
$\lm_3 = \lm_2$ and $\lm_5=\lm_4$, both in loop-momentum space and when
expanding the parametric integrals. So we are left with the NNLO Glauber
contributions yielding
\begin{multline}
  -i\pi^{d/2} \,
  \frac{i \, e^{-i\pi(\lm_1+2\lm_2+2\lm_4)} \,
      \Gm\left(\lm_2-\tfrac{1}{2}\right) \Gm\left(\lm_4-\tfrac{1}{2}\right)
      \Gm(\ep+\lm_1+2\lm_2+2\lm_4)}{
    16 (Q^2)^3 \, (m^2)^{\ep+\lm_1+2\lm_2+2\lm_4} \,
      \Gm(1+\lm_1) \Gm(1+\lm_2) \Gm(1+\lm_4) \Gm(-1-\ep)}
\\ \times
  \Gm(1-\ep-2\lm_2-2\lm_4)
  \label{Fp_l_NNLO_g}
\end{multline}
and the LO hard contribution,
\begin{multline}
  i\pi^{d/2} \,
  \frac{i \, e^{-i\pi(\lm_1+2\lm_2+2\lm_4)} \,
      e^{i\pi(\ep+\lm_1+2\lm_2+2\lm_4)/2} \,
      \Gm\left(\tfrac{-1-\ep-\lm_1-2\lm_2}{2}\right)
      \Gm\left(\tfrac{-1-\ep-\lm_1-2\lm_4}{2}\right)}{
    2\sqrt{\pi} \, (2Q^2)^{3+\ep+\lm_1+2\lm_2+2\lm_4} \,
      \Gm(1+\lm_2) \Gm(1+\lm_4) \Gm(-1-2\ep-\lm_1-2\lm_2-2\lm_4)}
\\ \times
  \Gm\left(\tfrac{-2-\ep-\lm_1-2\lm_2-2\lm_4}{2}\right)
  \Gm\left(\tfrac{3+\ep+\lm_1+2\lm_2+2\lm_4}{2}\right)
  \,,
  \label{Fp_l_NNLO_h}
\end{multline}
both consistent between the expansions of the loop integral and of the
parametric integrals. Setting all $\lm_i = 0$ in (\ref{Fp_l_NNLO_g})
and~(\ref{Fp_l_NNLO_h}), we recover the results from (\ref{Fp_NNLO_g})
and (\ref{Fp_NNLO_h}), respectively.

We may also evaluate the parametric integral~(\ref{Fp_integral_l_A}) in
terms of a onefold Mellin--Barnes representation:
\begin{align}
  I_{\pm} &= -i\pi^{d/2} \,
    \frac{e^{-i\pi(\lm_1+2\lm_2+2\lm_4)}}{
      4\pi \, \Gm(1+\lm_1) \Gm(1+\lm_2) \Gm(1+\lm_4)}
\nn \\  &\times
    \frac{1}{2\pi i} \int \dd z \, (m^2)^z
    (\pm 4 Q^2-i0)^{-3-\ep-\lm_1-2\lm_2-2\lm_4-z}
\nn \\  &\times
    \Gm(1+\lm_1+z) \Gm(3+\ep+\lm_1+2\lm_2+2\lm_4+z)
\nn \\  &\times
    \frac{\Gm(-z)
        \Gm\left(\tfrac{-1-\ep-\lm_1-2\lm_2-z}{2}\right)
        \Gm\left(\tfrac{-1-\ep-\lm_1-2\lm_4-z}{2}\right)
        \Gm^2\left(\tfrac{-2-\ep-\lm_1-2\lm_2-2\lm_4-z}{2}\right)}{
      \Gm(-1-2\ep-\lm_1-2\lm_2-2\lm_4-z)}
\;.
\label{Fp_integral_l_MB}
\end{align}
The relevant regions can easily be determined from the gamma functions in
the numerator of the last line. In particular, the squared gamma function
indicates that the expansion of the parametric integrals~$I_\pm$ requires
\emph{two} regions for the Glauber contribution, both scaling as
$(m^2)^{-2-\ep-\lm_1-2\lm_2-2\lm_4}$ at leading order.
%
When combining $F(Q^2,m^2) = 2(I_+ + I_-)$, one of these gamma functions is
cancelled, and we obtain
\begin{align}
  F(Q^2,m^2) &= i\pi^{d/2} \,
    \frac{i \, e^{-i\pi(\lm_1+2\lm_2+2\lm_4)}}{
      2\sqrt{\pi} \, \Gm(1+\lm_1) \Gm(1+\lm_2) \Gm(1+\lm_4)}
\nn \\  &\times
    \frac{1}{2\pi i} \int \dd z \, (m^2)^z
    (2 Q^2)^{-3-\ep-\lm_1-2\lm_2-2\lm_4-z} \,
    e^{i\pi(\ep+\lm_1+2\lm_2+2\lm_4+z)/2}
\nn \\  &\times
    \Gm(1+\lm_1+z) \Gm\left(\tfrac{3+\ep+\lm_1+2\lm_2+2\lm_4+z}{2}\right)
\nn \\  &\times
    \frac{\Gm(-z)
        \Gm\left(\tfrac{-1-\ep-\lm_1-2\lm_2-z}{2}\right)
        \Gm\left(\tfrac{-1-\ep-\lm_1-2\lm_4-z}{2}\right)
        \Gm\left(\tfrac{-2-\ep-\lm_1-2\lm_2-2\lm_4-z}{2}\right)}{
      \Gm(-1-2\ep-\lm_1-2\lm_2-2\lm_4-z)}
\;.
\label{Fp_integral_l_MB2}
\end{align}
%
{}From this representation the contributions to the asymptotic expansion in
the limit $m^2/Q^2 \to 0$ can be extracted: The hard contributions stem
from the residues of the poles at $z=n$, the 1-~and 2-collinear
contributions from $z=-1+2n-\ep-\lm_1-2\lm_2$ and
$z=-1+2n-\ep-\lm_1-2\lm_4$, respectively, the Glauber contributions from
$z=-2+2n-\ep-\lm_1-2\lm_2-2\lm_4$ (with $n=0,1,2,\ldots$). All LO, NLO and
NNLO contributions reported in (\ref{Fp_l_LO})--(\ref{Fp_l_NNLO_h}) are
confirmed by the corresponding residue contributions
from~(\ref{Fp_integral_l_MB2}). In particular, the structure of the poles
in~(\ref{Fp_integral_l_MB2}) clearly shows that the Glauber region does not
contribute to the NLO result, and that the collinear contributions are
absent at NNLO, as obtained before.

The results reported in this section show that the use of generic
propagator powers helps disentangling the Glauber and collinear regions
from each other, making all regions contribute in the same way to the
asymptotic expansion in loop-momentum space and to the expansion of the
parametric integrals. Keeping the dependence of the contributions and their
scalings on the propagator powers also facilitates the identification of
regions found by \texttt{asy2.m} in parametric space for a subsequent
expansion at the level of the loop integration.



\section{Summary of \texttt{asy2.m}}
\label{sec:syntax}

The updated version of the code, \texttt{asy2.m}, can be downloaded from
the known web site \cite{asy_web}, where further installation instructions
are found.
%
The \texttt{Mathematica} code is loaded using \texttt{<<asy2.m}.

The main function \texttt{AlphaRepExpand[]} reveals regions for loop
integrals:
\begin{verbatim}
  AlphaRepExpand[{k1, k2, ...},
    {(k1 + p1)^2 - m1^2, (k2 + p2)^2 - m2^2, ...},
    {p1^2 -> Q1, p2^2 -> Q2, p1*p2 -> Q3, ...},
    {m1^2 -> x, m2^2 -> x^2, Q1 -> 1, Q2 -> 3/2, ...},
    options]
\end{verbatim}
The first argument is the list of loop momenta. The second argument lists
the denominators of the propagators. The third argument contains
replacement rules for all kinematic invariants. In particular, all external
momenta appearing in the denominators must be replaced here, otherwise they
are not identified correctly as vectors. The fourth argument sets the
scaling of the parameters by replacing all symbols with powers of the
expansion parameter, labelled by the global symbol~\texttt{x}, and rational
numbers (integers or explicit fractions of integers).

The output is a list of regions, specified by the scaling (in powers of the
expansion parameter) of the Feynman parameters $x_1,x_2,\ldots$
corresponding to the propagators in the order stated in the second
argument.
%
Possible options are:
\begin{itemize}
\item \texttt{PreResolve -> True}: Try to eliminate cancellations between
  terms in the parametric representation by decomposing the integral and
  performing changes of variables. The code warns if this fails. The output
  contains, for each region, the changes of variables, the Jacobian of the
  integral transformation and the scalings of the new variables. Without
  this option, regions occurring at such cancellations will not be
  revealed.
\item \texttt{GenericPowers -> True}: Take into account generic (in
  particular non-integer) powers of the propagators, e.g.\ when these
  powers are used as analytic regulators. Without this option, the
  preresolution algorithm triggered by the option \texttt{PreResolve} only
  finds all regions if the propagator powers are positive integers.
\item \texttt{Verbose -> True}: Print verbose internal information.
\item \texttt{Scalar -> True}: Permit more complex structures of the
  denominators by specifying scalar products of momenta via the function
  \texttt{Scalar[k,p]} instead of simple products \texttt{k*p} or
  \texttt{Scalar2[k]} instead of \verb@k^2@, e.g.:
  \begin{verbatim}
AlphaRepExpand[{k},
  {Scalar[k, k] - m^2, Scalar2[k - q] - m^2},
  {Scalar2[q] -> qq, m^2 -> qq/4 + y}, {qq -> 1, y -> x}, 
  PreResolve -> True, Scalar -> True]
  \end{verbatim}
\end{itemize}

For expanding more general integrals, which need not originate from Feynman
diagrams, the command \texttt{WilsonExpand[]} may be used:
\begin{verbatim}
  WilsonExpand[F, U, {x1, x2, ...}, {... -> x, ...}, options]
\end{verbatim}
Traditionally, the first two arguments are the polynomials $\mathcal{F}$
and $\mathcal{U}$, respectively, from the parametric representation of the
Feynman integral. But, more generally, \texttt{WilsonExpand[]} will reveal
regions for integrals over parameters $x_1,x_2,\ldots$, integrated from 0
to~$\infty$ each, where all non-trivial polynomials of the parameters~$x_i$
occurring in the integrand are specified either in the first or the second
argument. The third argument of \texttt{WilsonExpand[]} is the list of
integration parameters, and the fourth argument specifies the scaling of
all quantities with the expansion parameter~\texttt{x}.

The output is a list of regions, specified by the scaling of the
integration parameters in powers of the expansion parameters.
%
Possible options are:
\begin{itemize}
\item \texttt{Delta -> True}: Under the integral, the sum over an arbitrary
  non-empty subset of the integration parameters is restricted to~1 via a
  delta function. The specific choice of the sum in this delta function
  must be irrelevant for the integral, which is the case for the
  generalized Feynman representations used in this paper.
\item \texttt{Verbose -> True}: Print verbose internal information.
\end{itemize}
More generally, the option \texttt{Delta} of \texttt{WilsonExpand[]} works
correctly for all integrals of the form
\begin{align}
  \int_0^\infty \cdots \int_0^\infty
  \dd x_1\cdots\dd x_n \,
  \delta\left(\sum_{i=1}^n a_i x_i - 1\right)
  f(x_1,\ldots,x_n)
  \;,
  \label{delta_int}
\end{align}
where the linear combination in the argument of the delta function has no
negative and at least one positive coefficient ($a_i \ge 0 \, \forall i$
and $\exists \, a_i>0$), and where the function~$f$ scales homogeneously
with the parameters~$x_i$ as
\begin{align}
  f(\lm x_1,\ldots,\lm x_n) = \lm^{-n} f(x_1,\ldots,x_n)
  \; \forall \lm>0
  \;,
  \label{delta_int_hom}
\end{align}
the degree of homogeneity being equal to minus the number of integration
parameters. It can be shown that for such integrals~(\ref{delta_int}) the
specific choice of the coefficients~$a_i$ is irrelevant.\footnote{%
  To see this, multiply the integrand of~(\ref{delta_int}) by~1 in the form
  $x_j \int_0^\infty \dd t \, e^{-t x_j}$, where $x_j$ is any of the
  integration parameters. Then, inside the $t$-integration, transform the
  other integration variables as $x_i \to x_i/t \, \forall i$ and use the
  homogeneity relation~(\ref{delta_int_hom}) with $\lm=1/t$. Finally
  evaluate the $t$-integration first, yielding $\int_0^\infty \dd t \,
  \delta(\sum_{i=1}^n a_i x_i - t) = 1$, independent of the
  coefficients~$a_i$, as long as the linear combination is positive. The
  integral~(\ref{delta_int}) is given by $\int_0^\infty \cdots
  \int_0^\infty \dd x_1\cdots\dd x_n \, f(x_1,\ldots,x_n) \, x_j \,
  e^{-x_j}$.}
%
All integrands of Feynman-parameter representations used in this paper
(without the additional analytic regularization factors $x_i^{\dl_i}$)
fulfill the homogeneity condition~(\ref{delta_int_hom}).



\section{Conclusion}
\label{sec:conclusion}

We have presented an algorithm for finding potential regions and Glauber
regions when expanding loop integrals at the level of their parametric
representations. The necessary decompositions and variable transformations
of the integral are automated by the updated \texttt{Mathematica} code
\texttt{asy2.m}. When the command \texttt{AlphaRepExpand[]} is used with
its option \texttt{PreResolve} enabled, we are convinced that it will
either reveal all relevant regions or issue a warning. In particular,
regions corresponding to cancellations between large positive and negative
terms in the parametric representation of the loop integrals (such as
potential and Glauber regions) will now be found.

Let us emphasize that to prove expansion by regions
at least for some specific limit typical of Minkowski
space is a natural mathematical problem. Perhaps, this problem is not specifically
related to Feynman integrals. Let us present an example of a one-dimensional
parametric integral, without any relevance to Feynman integrals, and show that
expansion by regions works successfully. To do this, we will use
{\tt asy2.m}.

Let us consider the integral
\bea
F(t)=\int_0^\infty(t + u + u^2)^\lm \dd u\;,
\label{toy_example}
\eea
with $\lm$ a complex parameter, in the limit $t\to 0$.
We assume that $\lm$ is in the domain Re$\lm<-1/2$ in order to have
an absolute convergence of the integral which then can be continued analytically
to the whole complex plane as an analytic function of $\lm$.
%
Running
\begin{verbatim}
  WilsonExpand[t + u + u^2, 1, {u}, {t -> x}]
\end{verbatim}
we obtain the two regions \verb@{{1}, {0}}@.
The leading-order terms from each region
can be evaluated analytically in terms of gamma functions at general $\lm$, with
the results
\be
\frac{t^{\lambda +1} \Gamma (-\lambda -1)}{\Gamma (-\lambda )}
\ee
and
\be
\frac{\Gamma (-2 \lambda -1) \Gamma (\lambda +1)}{\Gamma (-\lambda )} \;.
\ee
They can be checked easily by deriving the onefold Mellin--Barnes
representation
\be
F(t)=\frac{1}{2\pi i}\frac{1}{\Gamma (-\lambda )}\int  \Gamma (-z) \Gamma (\lambda -z+1) \Gamma (-2 \lambda +2
   z-1)t^z\dd z
\ee
and evaluating the first terms of the asymptotic expansion in the limit
$t\to 0$ by shifting the contour to the right and taking residues at the poles of
the two gamma functions in the integrand.



\section*{Acknowledgments}

This work is supported by the DFG Sonder\-forschungs\-bereich/Transregio~9
``Com\-pu\-ter\-gest\"utzte Theoretische Teilchenphysik''.
The work of A.S. and V.S. is also supported
by the Russian Foundation for Basic Research through grant 11-02-01196.
%
The authors thank M.~Beneke and A.~Pak for helpful discussions.



\begin{thebibliography}{99}

\bibitem{Chetyrkin:1988zz}
  K.G.~Chetyrkin,
  %``Operator Expansions In The Minimal Subtraction Scheme. 1: The Gluing Method,''
  Theor.\ Math.\ Phys.\  {\bf 75} (1988) 346
   [Teor.\ Mat.\ Fiz.\  {\bf 75} (1988) 26].
  %%CITATION = TMPHA,75,346;%%

\bibitem{Chetyrkin:1988cu}
  K.G.~Chetyrkin,
  %``Operator Expansions In The Minimal Subtraction Scheme. 2: Explicit Formulas For Coefficient Functions,''
  Theor.\ Math.\ Phys.\  {\bf 76} (1988) 809
   [Teor.\ Mat.\ Fiz.\  {\bf 76} (1988) 207].
  %%CITATION = TMPHA,76,809;%%

\bibitem{Gorishnii:1989dd}
  S.G.~Gorishny, %S.G.~Gorishnii,
  %``Construction Of Operator Expansions And Effective Theories In The Ms Scheme,''
  Nucl.\ Phys.\ B {\bf 319} (1989) 633.
  %%CITATION = NUPHA,B319,633;%%

\bibitem{Smirnov:1990rz}
  V.A.~Smirnov,
  %``Asymptotic expansions in limits of large momenta and masses,''
  Commun.\ Math.\ Phys.\  {\bf 134} (1990) 109.
  %%CITATION = CMPHA,134,109;%%

\bibitem{Smirnov:1994tg}
  V.A.~Smirnov,
  %``Asymptotic expansions in momenta and masses and calculation of Feynman diagrams,''
  Mod.\ Phys.\ Lett.\ A {\bf 10} (1995) 1485
  [hep-th/9412063].
  %%CITATION = HEP-TH/9412063;%%

\bibitem{books1a}%{Smirnov:2002pj}
  V.A.~Smirnov,
  {\em Applied asymptotic expansions in momenta and masses},
  STMP {\bf 177}, Springer, Berlin, Heidelberg (2002).
  %%CITATION = STPHB,177,1;%%

\bibitem{Seidensticker:1999bb}
  T.~Seidensticker,
  %``Automatic application of successive asymptotic expansions of Feynman diagrams,''
  hep-ph/9905298.
  %%CITATION = HEP-PH/9905298;%%

\bibitem{Harlander:1997zb}
  R.~Harlander, T.~Seidensticker and M.~Steinhauser,
  %``Complete corrections of O(alpha alpha(s)) to the decay of the Z boson  into
  %bottom quarks,''
  Phys.\ Lett.\ B {\bf 426} (1998) 125
  [hep-ph/9712228].
  %%CITATION = HEP-PH/9712228;%%

\bibitem{Beneke:1997zp}
  M.~Beneke and V.A.~Smirnov,
  %``Asymptotic expansion of Feynman integrals near threshold,''
  Nucl.\ Phys.\ B {\bf 522} (1998) 321
  [hep-ph/9711391].
  %%CITATION = HEP-PH/9711391;%%

\bibitem{Smirnov:1998vk}
  V.A.~Smirnov and E.R.~Rakhmetov,
  %``The Strategy of regions for asymptotic expansion of two loop vertex Feynman diagrams,''
  Theor.\ Math.\ Phys.\  {\bf 120} (1999) 870
   [Teor.\ Mat.\ Fiz.\  {\bf 120} (1999) 64]
  [hep-ph/9812529].
  %%CITATION = HEP-PH/9812529;%%

\bibitem{Smirnov:1999bza}
  V.A.~Smirnov,
  %``Problems of the strategy of regions,''
  Phys.\ Lett.\ B {\bf 465} (1999) 226
  [hep-ph/9907471].
  %%CITATION = HEP-PH/9907471;%%

\bibitem{Jantzen:2011nz}
  B.~Jantzen,
  %``Foundation and generalization of the expansion by regions,''
  JHEP {\bf 12} (2011) 076
  [arXiv:1111.2589 [hep-ph]].
  %%CITATION = ARXIV:1111.2589;%%

\bibitem{FIESTA2}
  A.V.~Smirnov, V.A.~Smirnov and M.~Tentyukov,
  %``FIESTA 2: parallelizeable multiloop numerical calculations,''
  Comput.\ Phys.\ Commun.\  {\bf 182} (2011) 790
  [arXiv:0912.0158 [hep-ph]].
  %%CITATION = CPHCB,182,790;%%

\bibitem{FIESTA}
  A.V.~Smirnov and M.N.~Tentyukov,
  %``Feynman Integral Evaluation by a Sector decomposiTion Approach (FIESTA),''
  Comput.\ Phys.\ Commun.\  {\bf 180} (2009) 735
  [arXiv:0807.4129 [hep-ph]].
  %%CITATION = CPHCB,180,735;%%

\bibitem{asy}%{Pak:2010pt}
  A.~Pak and A.~Smirnov,
  %``Geometric approach to asymptotic expansion of Feynman integrals,''
  Eur.\ Phys.\ J.\ C {\bf 71} (2011) 1626
  [arXiv:1011.4863 [hep-ph]].
  %%CITATION = ARXIV:1011.4863;%%

\bibitem{asy_web}
  \href{http://www-ttp.particle.uni-karlsruhe.de/~asmirnov/Tools-Regions.htm}
    {http://www-ttp.particle.uni-karlsruhe.de/\texttt{\symbol{126}}asmirnov/Tools-Regions.htm}

\end{thebibliography}

\end{document}

