path: root/doc/hhtop-doc/src/hhtop-doc.tex

\documentclass{kiss}
% load packages
\usepackage{header}
% BBlog bibliography commands
\usepackage{BBlog}
% miscellaneous commands
\usepackage{toolbox}
% main style file
\usepackage{iansecs}

\begin{document}
\hfil{\bf\Large hhtop}\par
\bigskip
\hfil{\bf v1.0}\par
\hugeskip

\indent {\tt hhtop} is a tool to compute, numerically, the following quantities for the Haldane-Hubbard model:
\begin{itemize}
\item the one-loop renormalization of the topological phase diagram,
\item the difference of the $(a,a)$ and $(b,b)$ wave-function renormalizations, at second order,
\end{itemize}
\hugeskip

\tableofcontents
\vfill\eject

\setcounter{page}1
\pagestyle{plain}

\section{Phase diagram}
\indent In this section we discuss the computation of the renormalization of the phase diagram.
\subseqskip

\subsection{Description of the computation}
\subsubsection{Definition of the problem}
\indent We wish to solve the following equation:
\begin{equation}
\tilde M_{\omega,t_1,t_2,\lambda}(W,\phi):=W+3\sqrt3\omega t_2\sin\phi+\frac{3\sqrt3}{16\pi^3}\lambda\int_{\mathcal B}dk\int_{-\infty}^\infty dk_0\ \frac{m_{t_2,W,\phi}(k)}{D_{t_1,t_2,W,\phi}(k_0,k)}=0
\label{eqrenmass}\end{equation}
for $W\in\mathbb{R}$, $\phi\in(-\pi,\pi]$, where the parameters $\omega=\pm1$, $t_2\geqslant0$, $t_1\geqslant3t_2$ and $\lambda\in\mathbb{R}$ are fixed. We now define the quantities appearing in~\-(\ref{eqrenmass}):
\begin{equation}
\mathcal B:=\left\{\left(\frac{2\pi}3+k'_1,k_2\right)\in\mathbb{R}^2\quad|\quad |k_2|<\frac{2\pi}{\sqrt3}-\sqrt3|k'_1|\right\},
\label{eqbrillouin}\end{equation}
\begin{equation}
\alpha_1(k_1,k_2):=\frac32+
\cos(\sqrt3k_2)+2\cos\left(\frac32k_1\right)\cos\left(\frac{\sqrt3}2k_2\right),
\label{eqalpha1}\end{equation}

\begin{equation}
\alpha_2(k_1,k_2):=
-\sin(\sqrt3k_2)+2\cos\left(\frac32k_1\right)\sin\left(\frac{\sqrt3}2k_2\right),
\label{eqalpha2}\end{equation}

\begin{equation}
\Omega(k_1,k_2):=1+2e^{-\frac32ik_1}\cos\left(\frac{\sqrt3}2k_2\right),
\label{eqOmega}\end{equation}

\begin{equation}
m_{t_2,W,\phi}(k):=W-2t_2\sin\phi\alpha_2(k)
\label{eqm}\end{equation}

\begin{equation}
\zeta_{t_2,\phi}(k):=2t_2\cos\phi\alpha_1(k),\quad
\xi_{t_1,t_2,W,\phi}(k):=\sqrt{m_{t_2,W,\phi}^2(k)+t_1^2|\Omega(k)|^2}
\label{eqOmega}\end{equation}

\begin{equation}
D_{t_1,t_2,W,\phi}(k_0,k):=(ik_0+\zeta_{t_2,\phi}(k))^2-\xi_{t_1,t_2,W,\phi}^2(k).
\label{eqDdef}\end{equation}
\bigskip

\subsubsection{Integration of the Matsubara momentum}
\indent We first integrate out $k_0$ analytically. We use the following identity: for $x\in\mathbb{R}$ and $y>0$,
\begin{equation}
\int_{-\infty}^\infty dk_0\ \frac1{((ik_0+x)^2-y^2}=-\chi(x^2<y^2)\frac\pi y
\label{eqintk0}\end{equation}
in which $\chi(x^2<y^2)\in\{1,0\}$ is equal to 1 if and only if $x^2<y^2$. Furthermore (see appendix~\-(\ref{appintk0})), if
\begin{equation}
t_1\geqslant3t_2
\label{eqcondt}\end{equation}
then
\begin{equation}
\zeta_{t_2,\phi}^2(k)\leqslant\xi_{t_1,t_2,W,\phi}^2(k)
\label{eqineqab}\end{equation}
for all $k\in\mathcal B$, $\phi\in(-\pi,\pi]$ and $W\in\mathbb{R}$, which implies that
\begin{equation}
\tilde M_{\omega,t_1,t_2,\lambda}(W,\phi)=W+3\sqrt3\omega t_2\sin\phi-\frac{3\sqrt3}{16\pi^2}\lambda\int_{\mathcal B}dk\ \frac{m_{t_2,W,\phi}(k)}{\xi_{t_1,t_2,W,\phi}(k)}.
\label{eqrenmassintk0}\end{equation}
\bigskip

\subsubsection{Reduction by symmetries}
\indent By using some symmetries of the integrand of~\-(\ref{eqrenmassintk0}), we can reduce the integration region. Indeed, $m_{t_2,W,\phi}(k)$ and $\xi_{t_1,t_2,W,\phi}(k)$ are symmetric under $k_1\mapsto-k_1$ and under rotations of angle $\frac{2\pi}3$. In addition, $m_{t_2,W,\phi}(k_1,k_2)=m_{t_2,W,-\phi}(k_1,-k_2)$ and $\xi_{t_1,t_2,W,\phi}(k_1,k_2)=\xi_{t_1,t_2,W,-\phi}(k_1,-k_2)$. We can therefore rewrite
\begin{equation}
\tilde M_{\omega, t_1, t_2,\lambda}(W,\phi)=W+3\sqrt3\omega t_2\sin\phi-\lambda(I_{t_1,t_2}(W,\phi)+I_{t_1,t_2}(W,-\phi))
\label{eqrenmassintk0half}\end{equation}
where
\begin{equation}
I_{t_1,t_2}(W,\phi):=
\frac{9\sqrt3}{8\pi^2}\int_{\mathcal B_+}dk\ \frac{m_{t_2,W,\phi}(k)}{\xi_{t_1,t_2,W,\phi}(k)}.
\label{eqI}\end{equation}
with 
\begin{equation}
\mathcal B_+:=\left\{(k_1,k_2)\in\mathcal B\ |\ k_2>0,\ k_1<\frac23\pi,\ k_2<\frac1{\sqrt3}k_1\right\}.
\label{eqB0}\end{equation}
\bigskip

\subsubsection{Polar coordinates}
\indent Let
\begin{equation}
p_F^\pm:=\left(\frac{2\pi}3,\ \frac{2\pi}{3\sqrt3}\right).
\label{eqfermi}\end{equation}
We note that $\xi_{t_1,t_2,W,\phi}$ has roots if and only if $m_{t_2,W,\phi}(p_F^+)=0$ or $m_{t_2,W,\phi}(p_F^-)=0$, located at $p_F^+$ in the former case and at $p_F^-$ in the latter. If $m_{t_2,W,\phi}$ vanishes at both $p_F^\pm$, which can only occur if $W=0$ and $\phi=0,\pi$, then $\xi_{t_1,t_2,W,\phi}$ vanishes at both $p_F^\pm$. Nevertheless, the integrand in~\-(\ref{eqI}) is not singular, since $\xi_{t_1,t_2,W,\phi}(k'+p_F^+)\sim t_1|k'|$, and the integration over $k$ is 2-dimensional. In order to make this lack of singularity apparent, it is convenient to switch to polar coordinates around $p_F^+$: $(k_1,k_2)=p_F^++\frac{2\pi}{3\sqrt3}\rho(\cos\theta,\sin\theta)$:
\begin{equation}
I_{t_1,t_2}(W,\phi)=\frac{\sqrt3}{6}\int_{-\frac\pi6}^{\frac\pi6} d\theta\int_0^{R(\theta)} d\rho\ \rho\frac{\bar m_{t_2,W,\phi}(\rho,\theta)}{\bar\xi_{t_1,t_2,W,\phi}(\rho,\theta)},
\label{eqI}\end{equation}
in which
\begin{equation}\begin{array}c
\bar m_{t_2,W,\phi}(\rho,\theta):=W-2t_2\sin\phi\bar\alpha_2(\rho,\theta),\quad
\bar\xi_{t_1,t_2,W,\phi}(\rho,\theta):=\sqrt{m_{t_2,W,\phi}^2(\rho,\theta)+t_1^2|\bar\Omega(\rho,\theta)|^2}\\[0.5cm]
\bar\alpha_2(\rho,\theta):=-2\sin\left(\frac\pi3(1+\rho\sin\theta)\right)\left(\cos\left(\frac\pi3(1+\rho\sin\theta)\right)+\cos\left(\frac\pi{\sqrt3}\rho\cos\theta\right)\right)\\[0.5cm]
|\bar\Omega(\rho,\theta)|^2=1+4\cos\left(\frac\pi3(1+\rho\sin\theta)\right)\left(\cos\left(\frac\pi3(1+\rho\sin\theta)\right)-\cos\left(\frac\pi{\sqrt3}\rho\cos\theta\right)\right)
\end{array}\label{eqxipolar}\end{equation}
and
\begin{equation}
R(\theta):=\frac1{\cos(\theta-\frac\pi6)}.
\label{eqpolarbound}\end{equation}

\vskip10pt
\subsection{Strategy of the numerical computation}
\subsubsection{Newton scheme}\label{secnewtonphase}
\indent In order to solve~\-(\ref{eqrenmass}), we will use a Newton scheme (see section~\-\ref{secnewton}). More precisely, we fix $\phi$ and compute $W(\phi)$ as the limit of
\begin{equation}\begin{array}{>\displaystyle c}
W_0(\phi):=-\omega 3\sqrt3t_2\sin\phi\\[0.3cm]
W_{n+1}(\phi):=W_n(\phi)-\frac{\tilde M_{\omega,t_1,t_2,\lambda}(W_n(\phi),\phi)}{\partial_W\tilde M_{\omega,t_1,t_2,\lambda}(W_n(\phi),\phi)}.
\end{array}\label{eqnewton}\end{equation}
The first two derivatives of $\tilde M$ are
\begin{equation}
\partial_W\tilde M_{\omega,t_1,t_2,\lambda}(W,\phi)=1-\lambda(\partial_WI_{t_1,t_2}(W,\phi)+\partial_WI_{t_1,t_2}(W,-\phi))
\label{eqdM}\end{equation}
with
\begin{equation}
\partial_WI_{t_1,t_2}(W,\phi)=\frac{\sqrt3}{6}\int_{-\frac\pi6}^{\frac\pi6}d\theta\int_0^{R(\theta)}d\rho\ \frac{\rho}{\bar\xi_{t_1,t_2,W,\phi}(\rho,\theta)}\left(1-\frac{\bar m^2_{t_2,W,\phi}(\rho,\theta)}{\bar\xi_{t_1,t_2,W,\phi}^2(\rho,\theta)}\right)
\label{eqdI}\end{equation}
and
\begin{equation}
\partial_W^2\tilde M_{\omega,t_1,t_2,\lambda}(W,\phi)=-\lambda(\partial_W^2I_{t_1,t_2}(W,\phi)+\partial_W^2I_{t_1,t_2}(W,-\phi))
\label{eqddM}\end{equation}
with
\begin{equation}
\partial_W^2I_{t_1,t_2}(W,\phi)=-\frac{\sqrt3}{2}\int_{-\frac\pi6}^{\frac\pi6}d\theta\int_0^{R(\theta)}d\rho\ \frac{\rho\bar m_{t_2,W,\phi}(\rho,\theta)}{\bar\xi^3_{t_1,t_2,W,\phi}(\rho,\theta)}\left(1-\frac{\bar m^2_{t_2,W,\phi}(\rho,\theta)}{\bar\xi_{t_1,t_2,W,\phi}^2(\rho,\theta)}\right).
\label{eqddI}\end{equation}

\subsubsection{Integration}\label{secintegrationphase}
\indent In order to compute $W_n(\phi)$, we have to evaluate $I_{t_1,t_2}(W,\phi)$ and $\partial_W I_{t_1,t_2}(W,\phi)$. The integrations are carried out using Gauss-Legendre quadratures (see section~\-\ref{secintegrationgl}). In order to use this method to compute the double integral over $\theta$ and $\rho$, we rewrite
\begin{equation}
\int d\theta\int d\rho\ F(\theta,\rho)=\int d\theta\ G(\theta),\quad
G(\theta):=\int d\rho\ F(\theta,\rho)
\label{eqmultidimint}\end{equation}
for the appropriate $F$.

\subsection{Usage and examples}
\indent We will now describe some basic usage cases of {\tt hhtop phase}. For a full description of the options of {\tt hhtop}, see the {\tt man} page.
\subseqskip

\subsubsection{Basic usage}
\indent The value of the parameters can be set via the {\tt -p} flag. Here is an example\par
\medskip
\indent{\tt hhtop phase -p "omega=1;t1=1.;t2=.1;lambda=.01;sinphi=1;"}\par
\medskip
Note that $\phi$ can be set instead of $\sin\phi$, though the result of the computation only depends on $\sin\phi$. The parameters that are not specified by the {\tt-p} flag are set to their default value: $\omega=1$, $t_1=1$, $t_2=0.1$, $\lambda=0.01$, $\sin\phi=1$.

\subsubsection{Precision of the computation}
\indent The precision of the computation can be controlled by three parameters: the precision of the numbers manipulated by {\tt hhtop} (set via the {\tt -P} flag, see section~\-\ref{secprecision}), the order of the integration, and the tolerance of the Newton scheme.
\bigskip

\point{\bf Order of the integration.} The order of the integration, that is, the value of the number $N$ introduced in section~\-\ref{secintegrationgl}, can be specified via the {\tt -O} flag. Its default value is 10. The difference of the value of the integral at different orders is a good measure of the numerical error. Example:\par
\medskip
\indent{\tt hhtop phase -O 30}
\bigskip

\point{\bf Tolerance of the Newton scheme.} The Newton iteration halts when the difference $|x_{n+1}-x_n|$ (see section~\-\ref{secnewton}) is smaller than a number, called the {\it tolerance} of the algorithm, or when the toal number of steps exceeds a given threshold. The tolerance can be set via the {\tt-t} flag, and the maximal number of steps via the {\tt -N} flag. Their default values are $10^{-11}$ and $1000000$. The tolerance and maximal number of steps are also used in the computation of the roots $\{x_1,\cdots,x_N\}$ of the $N$-th Legendre polynomial which are used for the numerical integration (see section~\ref{secintegrationgl}). If the tolerance or the maximal number of steps are too small, and the precision of multi-precision floats is too low, then the iteration may not converge. Example:\par
\medskip
\indent{\tt hhtop phase -t 1e-30 -N 2000000 -O 100 -P 256}

\subsubsection{Using double precision floats instead of multi-precision floats}
\indent Using the {\tt -D} command-line flag, {\tt hhtop} can be instructed to use {\tt long double}'s instead of MPFR floats. Whereas one then loses the ability of adjusting the precision, the computation time can be drastically reduced. Example:\par
\medskip
\indent{\tt hhtop -D phase -p "sinphi=1.;"}\par
\bigskip
The precision of {\tt long double}s is compiler dependent, see section~\-\ref{secprecision}.

\vfill\eject

\section{Wave function renormalization}\label{secz1z2}
\indent In this section we discuss the computation of the difference and the sum of the $(a,a)$ and $(b,b)$ wave-function renormalizations.
\bigskip

{\bf Warning: This computation is only accurate if $\phi$ is not too close to $0$.}
\subseqskip

\subsection{Description of the computation}
\subsubsection{Definition of the problem}
\indent We wish to compute the following quantities:
\begin{equation}
  \begin{array}{>\displaystyle c}
    z_1-z_2=i\frac{27}{128\pi^4}(\partial_{k_0}S_+|_{k_0=0}-\partial_{k_0}S_-|_{k_0=0})
    \\[0.5cm]
    z_1+z_2=i\frac{27}{128\pi^4}(\partial_{k_0}S_+|_{k_0=0}+\partial_{k_0}S_-|_{k_0=0})
  \end{array}
  \label{eqz1z2}
\end{equation}
where
\begin{equation}
S_\pm(k_0)=\int_{\mathcal B} dpdq\ \int_{-\infty}^\infty \frac{dp_0 dq_0}{2\pi^2}\ \frac{(-ip_0-\zeta_p\mp m_p)(-iq_0-\zeta_q\mp m_q)(-i(p_0+q_0-k_0)-\zeta_F\mp m_F)}{((ip_0+\zeta_p)^2-\xi_p^2)((iq_0+\zeta_q)^2-\xi_q^2)((i(p_0+q_0-k_0)+\zeta_F)^2-\xi_F^2)}
\label{eqS}\end{equation}
in which
\begin{equation}
\mathcal B:=\left\{\left(\frac{2\pi}3+k'_1,k_2\right)\in\mathbb{R}^2\quad|\quad |k_2|<\frac{2\pi}{\sqrt3}-\sqrt3|k'_1|\right\},
\label{eqbrillouin}\end{equation}
\begin{equation}
\alpha_1(k_1,k_2):=
2\cos\left(\frac{\sqrt3}2k_2\right)\left(\cos\left(\frac32k_1\right)+\cos\left(\frac{\sqrt3}2k_2\right)\right)+\frac12,
\label{eqalpha1}\end{equation}
\begin{equation}
\alpha_2(k_1,k_2):=
2\sin\left(\frac{\sqrt3}2k_2\right)\left(\cos\left(\frac32k_1\right)-\cos\left(\frac{\sqrt3}2k_2\right)\right),
\label{eqalpha2}\end{equation}
\begin{equation}
m(k):=W-2t_2\sin\phi\alpha_2(k)
\label{eqm}\end{equation}
\begin{equation}
\zeta(k):=2t_2\cos\phi\alpha_1(k),\quad
\xi(k):=\sqrt{m^2(k)+2t_1^2\alpha_1(k)}
\label{eqOmega}\end{equation}
and
\begin{equation}\begin{array}c
\zeta_p\equiv\zeta(p),\quad
\zeta_q\equiv\zeta(q),\quad
\zeta_F\equiv\zeta(p+q-p_F^\omega),\quad
\xi_p\equiv\xi(p),\quad
\xi_q\equiv\xi(q),\quad
\xi_F\equiv\xi(p+q-p_F^\omega),\\[0.3cm]
m_p\equiv m(p),\quad
m_q\equiv m(q),\quad
m_F\equiv m(p+q-p_F^\omega)
\end{array}\label{eqxip}\end{equation}
with $\omega\in\{-1,+1\}$ and
\begin{equation}
  p_F^\pm:=\left(\frac{2\pi}3,\pm\frac{2\pi}{3\sqrt3}\right).
\end{equation}

\subsubsection{Integration of the Matsubara momentum}
\indent We first integrate out $p_0$ and $q_0$ analytically.
We recall (see appendix~\-(\ref{appintk0})) that, provided $t_1\geqslant3t_2$,
\begin{equation}
\zeta^2(k)\leqslant\xi^2(k).
\label{eqineqab2}\end{equation}
By closing the integration path over $p_0$ around the positive-imaginary half-plane (which, by~\-(\ref{eqineqab2}), contains two poles), and using the residues theorem, we find that
\begin{equation}\begin{largearray}
S_\pm(k_0)=-\int_{\mathcal B} dpdq\int_{-\infty}^\infty \frac{dq_0}{2\pi}\ \left(\frac{(\xi_p\mp m_p)(-iq_0-\zeta_q\mp m_q)(-i(q_0-k_0)+\zeta_p-\zeta_F+\xi_p\mp m_F)}{\xi_p((iq_0+\zeta_q)^2-\xi_q^2)((i(q_0-k_0)-\zeta_p+\zeta_F-\xi_p)^2-\xi_F^2)}\right.\\[0.5cm]
\hfill+\left.\frac{(i(q_0-k_0)-\zeta_p+\zeta_F+\xi_F\mp m_p)(-iq_0-\zeta_q\mp m_q)(\xi_F\mp m_F)}{((-i(q_0-k_0)+\zeta_p-\zeta_F-\xi_F)^2-\xi_p^2)((iq_0+\zeta_q)^2-\xi_q^2)\xi_F}\right).
\end{largearray}\label{eqintmatsubara1}\end{equation}
We then close the integration path over $q_0$ around the positive-imaginary half-plane for the first term, and the negative imaginary half-plane for the second, and find
\begin{equation}\begin{array}{>\displaystyle r@{\ }>\displaystyle l}
S_\pm(k_0)=\frac12\int_{\mathcal B} dpdq&\left(
\frac{(\xi_p\mp m_p)(\xi_q\mp m_q)(ik_0+Z+\xi_p+\xi_q\mp m_F)}{\xi_p\xi_q((ik_0+Z+\xi_p+\xi_q)^2-\xi_F^2)}
\right.\\[0.5cm]&+\left.
\frac{(\xi_q\pm m_q)(\xi_F\mp m_F)(ik_0+Z-\xi_q-\xi_F\pm m_p)}{\xi_q\xi_F((ik_0+Z-\xi_q-\xi_F)^2-\xi_p^2)}
\right.\\[0.5cm]&-\left.
\frac{(\xi_p\mp m_p)(\xi_F\mp m_F)(ik_0+Z+\xi_p-\xi_F\pm m_q)}{\xi_p\xi_F((ik_0+Z+\xi_p-\xi_F)^2-\xi_q^2)}
\right)
\end{array}\label{eqSnok0}\end{equation}
with
\begin{equation}
Z:=\zeta_p+\zeta_q-\zeta_F.
\label{eqZ}\end{equation}
The sum of the three terms in the right side of~\-(\ref{eqSnok0}) yields
\begin{equation}\begin{largearray}
S_\pm(k_0)=\frac12\int_{\mathcal B} dpdq\ \left(\frac{(\xi_p\xi_q\xi_F-(\xi_pm_q+\xi_qm_p)m_F+m_pm_q\xi_F)(ik_0+Z)}{\xi_p\xi_q\xi_F((ik_0+Z)^2-(\xi_p+\xi_q+\xi_F)^2)}\right.\\[0.5cm]
\hfill\left.\pm\frac{(\xi_p+\xi_q+\xi_F)(m_p\xi_q\xi_F+\xi_pm_q\xi_F-\xi_p\xi_qm_F-m_pm_qm_F)}{\xi_p\xi_q\xi_F((ik_0+Z)^2-(\xi_p+\xi_q+\xi_F)^2)}\right).
\end{largearray}\label{eqSnok0simp}\end{equation}
Therefore,
\begin{equation}
  \begin{array}{>\displaystyle c}
    z_1-z_2
    =\frac{27}{64\pi^4}\int_{\mathcal B} dpdq\ \left(\frac{(\xi_p+\xi_q+\xi_F)(\frac{m_p}{\xi_p}+\frac{m_q}{\xi_q}-\frac{m_F}{\xi_F}-\frac{m_pm_qm_F}{\xi_p\xi_q\xi_F})Z}{(Z^2-(\xi_p+\xi_q+\xi_F)^2)^2}\right).
    \\[0.5cm]
    z_1+z_2
    =\frac{27}{128\pi^4}\int_{\mathcal B} dpdq\ \left(\frac{\left(1-\frac{m_pm_F}{\xi_p\xi_F}-\frac{m_qm_F}{\xi_q\xi_F}+\frac{m_pm_q}{\xi_p\xi_q}\right)(Z^2+(\xi_p+\xi_q+\xi_F)^2)}{(Z^2-(\xi_p+\xi_q+\xi_F)^2)^2}\right).
  \end{array}
  \label{eqz1z2}
\end{equation}

\subsubsection{Singularities of the integrand}
\indent In order to compute the integrals in~\-(\ref{eqz1z2}) numerically, we will use Gauss quadratures, which are only accurate if the integrands are smooth (i.e. if high order derivatives of the integrand are bounded). In this case, the integrand has singularities, indeed
\begin{itemize}
  \item $\alpha_1$ and $\zeta$ vanish at $p_F^+$ and $p_F^-$, and if $W=-\omega3\sqrt3t_2\sin\phi$, then $m$ vanishes at $p_F^\omega$,
  \item if $W=-\omega3\sqrt3t_2\sin\phi$, then the second derivative of $\xi$ diverges at $p_F^\omega$.
\end{itemize}
The asymptotics near the singularities are
\begin{equation}
  \begin{array}{r@{\ }>\displaystyle l}
    \sqrt{2t_1^2\alpha_1(p_F^\omega+k')}=&\frac32t_1|k'|+t_1\ O(|k'|^2)\\[0.3cm]
    \zeta(p_F^\omega+k')=&t_2\cos\phi\ O(|k'|^2)\\[0.3cm]
    m(p_F^\omega+k')-(W-\omega3\sqrt3t_2\sin\phi)=&t_2\sin\phi\ O(|k'|^2)
  \end{array}
\end{equation}
which implies that, if $W=-\omega3\sqrt3t_2\sin\phi$, $p=p_F^\omega+p'$, $q=p_F^\omega+q'$ and $k=p_F^\omega$, then
\begin{equation}
  \begin{array}{>\displaystyle c}
    \xi_p+\xi_q+\xi_F=\frac32t_1(|p'|+|q'|+|p'+q'|)(1+O(|p'|)+O(|q'|)),\\[0.3cm]
    Z=O(|p'|^2)+O(|q'|^2)+O(|p'+q'|^2),\quad
    \frac{m_p}{\xi_p}=O(|p'|),\quad
    \frac{m_q}{\xi_q}=O(|q'|),\quad
    \frac{m_F}{\xi_F}=O(|p'+q'|).
  \end{array}
  \label{eqasymp}
\end{equation}
In addition, the $O(\cdot)$ factors in~\-(\ref{eqasymp}) are analytic functions of $|p'|$, $|q'|$ and $|p'+q'|$. Note that, since $|\cdot|$ is not an analytic function (its second derivative diverges at $0$), the $O(\cdot)$ factors are {\it not} analytic functions of $p'$, $q'$ or $p'+q'$. Therefore, if $W=-\omega 3\sqrt3 t_2\sin\phi$, then
\begin{equation}
  \mathcal I_-(p,q):=\frac{27}{64\pi^4}\frac{(\xi_p+\xi_q+\xi_F)(\frac{m_p}{\xi_p}+\frac{m_q}{\xi_q}-\frac{m_F}{\xi_F}-\frac{m_pm_qm_F}{\xi_p\xi_q\xi_F})Z}{(Z^2-(\xi_p+\xi_q+\xi_F)^2)^2}
  \label{eqintegrand-}
\end{equation}
\begin{itemize}
  \item is smooth as long as $p\neq p_F^\omega$ and $q\neq p_F^\omega$ and $p+q\neq 2p_F^\omega$,
  \item is bounded for all $p,q$,
  \item its derivatives diverge if $p=p_F^\omega$ or $q=p_F^\omega$ or $p+q=2p_F^\omega$.
\end{itemize}
Similarly, if $W=-\omega 3\sqrt3 t_2\sin\phi$, then
\begin{equation}
  \mathcal I_+(p,q):=\frac{27}{128\pi^4}\frac{\left(1-\frac{m_pm_F}{\xi_p\xi_F}-\frac{m_qm_F}{\xi_q\xi_F}+\frac{m_pm_q}{\xi_p\xi_q}\right)(Z^2+(\xi_p+\xi_q+\xi_F)^2)}{(Z^2-(\xi_p+\xi_q+\xi_F)^2)^2}
  \label{eqintegrand+}
\end{equation}
\begin{itemize}
  \item is smooth as long as $p\neq p_F^\omega$ and $q\neq p_F^\omega$ and $p+q\neq 2p_F^\omega$,
  \item diverges if $p=p_F^\omega$ and $q=p_F^\omega$ (it would remain bounded if it were multiplied by $|p-p_F^\omega|\cdot|q-p_F^\omega|$),
  \item is bounded  for all $(p,q)\neq(p_F^\omega,p_F^\omega)$,
  \item its derivatives diverge if $p=p_F^\omega$ or $q=p_F^\omega$ or $p+q=2p_F^\omega$.
\end{itemize}
In the next section, we will regularize these singularities by changing performing an appropriate change of variables.

\subsubsection{Sunrise coordinates}
\indent In this section, we will show how to regularize the singularities mentioned in the previous section. We assume throughout this section that $W=-\omega 3\sqrt3t_2\sin\phi$ (if this is not the case, then there are no singularities).
\bigskip

{\bf Warning}: As it is set up here, {\bf this computation is only accurate if $\phi$ is not too close to 0} (see the remark on p.~\-\ref{rkphi}).
\bigskip

\indent While $\mathcal I_-$ and $|p-p_F^\omega||q-p_F^\omega|\mathcal I_+$ are singular functions of $p$ and $q$ (because of the divergence of the second derivative of $|p-p_F^\omega|$), they can be re-expressed as smooth functions of $p$, $q$, $\rho:=|p-p_F^\omega|$, $r:=|q-p_F^\omega|$ and $\gamma:=|p+q-2p_F^\omega|$. We will, therefore, change to the {\it sunrise} coordinates, described in appendix~\-\ref{appsunrise}, which, by lemma~\-\ref{lemmasunrise}, regularize the singularities of $\mathcal I_-$ and $\mathcal I_+$. However, the sunrise coordinates are only defined for rotationally symmetric integration regions, so we will have to split the integration regions, and note that in the regions where we cannot change to sunrise coordinates, it suffices to use polar coordinates.
\bigskip

\indent Let
\begin{equation}
  \mathcal B_\pm:=\mathcal B\cap\{(k_1,k_2)\in\mathcal B\ |\ \pm k_2>0\}
\end{equation}
and
\begin{equation}
  \mathcal B^{(F)}_\pm:=\left\{p_F^\pm+k',\ |k'|<R\right\},\quad
  R:=\frac{2\pi}{3\sqrt3},\quad
  \mathcal B^{(R)}_\pm:=\mathcal B_\pm\setminus\mathcal B^{(F)}_\pm
\end{equation}
($\mathcal B^{(F)}_\pm$ is the largest disk that is included in $\mathcal B_\pm$). As is discussed below (see the remark on p.~\-\ref{rkcutoff}), it is inconvenient to sharply split the integral, so we will use a smooth cut-off function instead: we define, for $\tau\in(0,1)$, $\chi_\tau:[0,\infty)\to[0,1]$:
\begin{equation}
  \chi_\tau(x):=
  \left\{\begin{array}{l@{\quad}>\displaystyle l}
    \frac{e^{-\frac{1-\tau}{1-x}}}{e^{-\frac{1-\tau}{x-\tau}}+e^{-\frac{1-\tau}{1-x}}(1-e^{-\frac{1-\tau}{x-\tau}})}&\mathrm{if\ }x\in(\tau,1)\\[0.5cm]
    0&\mathrm{if\ }x\in[0,\tau]\cup[1,\infty)
  \end{array}\right.
\end{equation}
which is equal to $1$ if $x\leqslant\tau$, and to $0$ if $x\geqslant1$, and is $\mathcal C^\infty$. In addition, one can prove that $\chi_\tau$ is a class-2 Gevrey function, that is, $\exists C_0,C>0$ such that for all $x\in[0,\infty)$ and $n\in\mathbb N$,
\begin{equation}
  \mathrm{sup}\left|\frac{d^n\chi_\tau}{dx^n}\right|\leqslant C_0C^n(n!)^2.
\end{equation}
Note that $C_0$ and $C$ depend on $\tau$, and diverge as $\tau\to1$. We will fix $\tau=\frac12$ in the following.
\bigskip

{\bf Remark}: By introducing such a cutoff function, the integrands will no longer be a analytic, but class-2 Gevrey functions. By lemma~\-\ref{lemmaGL} (see appendix~\-\ref{appGL}), the error of the numerical integration scheme nevertheless decays as an exponential in $\sqrt N$ where $N$ denotes the order of the quadrature.
\bigskip

Let, for $p\in\mathcal B$
\begin{equation}
  f_{\omega}^{(F)}(p):=\chi_{\frac12}\left(\frac{|p-p_F^\omega|_{\mathcal B}}R\right)\quad
  f_{\omega}^{(R)}(p):=1-f_{\omega}^{(F)}(p)
\end{equation}
where the choice $\tau=\frac12$ is arbitrary (any other value would do, as long as it is not too close to $0$ or $1$), we recall that $R:=\frac{2\pi}{3\sqrt3}$, and $|\cdot|_{\mathcal B}$ denotes the periodic Euclidian norm on $\mathcal B$:
\begin{equation}
  |k|_{\mathcal B}:=\min\{|k+n_1G_++n_2G_-|,\ (n_1,n_2)\in\mathbb Z^2\},\quad
  \textstyle G_\pm:=(\frac{2\pi}3,\pm\frac{2\pi}{\sqrt3}).
\end{equation}
We then split
\begin{equation}
    z_1\mp z_2=A_{F,F}^{(\mp)}+2A_{R,F}^{(\mp)}+A_{R,R}^{(\mp)}
\end{equation}
where
\begin{equation}
  \begin{array}{r@{\ }>\displaystyle l}
    A_{F,F}^{(\mp)}:=&\int_{\mathcal B_\omega^{(F)}}dp\int_{\mathcal B_\omega^{(F)}}dq\ f_{\omega}^{(F)}(p)f_{\omega}^{(F)}(q)\ \mathcal I_\mp(p,q)\\[0.5cm]
    A_{R,F}^{(\mp)}:=&\int_{\mathcal B}dp\int_{\mathcal B^{(F)}_\omega}dq\ f_{\omega}^{(R)}(p)f_{\omega}^{(F)}(q)\ \mathcal I_\mp(p,q)\\[0.5cm]
    A_{R,R}^{(\mp)}:=&\int_{\mathcal B}dp\int_{\mathcal B}dq\ f_{\omega}^{(R)}(p)f_{\omega}^{(R)}(q)\ \mathcal I_\mp(p,q)
  \end{array}
\end{equation}
in which we used the symmetry $\mathcal I_\mp(p,q)=\mathcal I_\mp(q,p)$. We will change to sunrise coordinates in $A_{F,F}$ and to polar coordinates in $A_{R,F}$ and $A_{R,R}$.
\bigskip

\point The integrand of $A_{F,F}^{(\mp)}$ has the same singularities as $\mathcal I_{\mp}$, which we regularize by changing to {\it sunrise coordinates}. Since $\mathcal B_\omega^{(F)}$ is a disk, these coordinates are well defined (see lemma~\-\ref{lemmasunrise}). In order to get rid of factors of $\pi$, we first rescale $p$ and $q$ by $R=\frac{2\pi}{3\sqrt3}$, and find
\begin{equation}
  A_{F,F}^{(\mp)}=2\int_0^1d\rho\int_0^{2\pi}d\theta\int_{-\frac\pi2}^{\frac\pi2}d\psi\int_0^1dz\ \Sigma f_{\omega,1}^{(F)}(\sigma)\Sigma f_{\omega,2}^{(F)}(\sigma)\Sigma J(\sigma)\Sigma\mathcal I_\mp(\sigma)
\end{equation}
with $\sigma\equiv(\rho,\theta,\psi,z)$,
\begin{equation}
  \Sigma J(\sigma)=
  4\rho^3\frac{\bar r(1+\bar r\cos(2\psi))}{(1+\cos\psi)\sqrt{1+\bar r\cos^2\psi}},
\end{equation}
where $\bar r$ is defined in~\-(\ref{eqr}),
\begin{equation}
  \Sigma f_{\omega,1}^{(F)}(\sigma):=\chi_{\frac12}\left(\rho\right),\quad
  \Sigma f_{\omega,2}^{(F)}(\sigma):=\chi_{\frac12}\left(\rho\bar r\right),
\end{equation}
\begin{equation}
  \begin{array}{>\displaystyle c}
    \Sigma\mathcal I_-(\sigma):=\frac{1}{108}\frac{(\Sigma\xi_p+\Sigma\xi_q+\Sigma\xi_F)(\frac{\Sigma m_p}{\Sigma\xi_p}+\frac{\Sigma m_q}{\Sigma\xi_q}-\frac{\Sigma m_F}{\Sigma\xi_F}-\frac{\Sigma m_p\Sigma m_q\Sigma m_F}{\Sigma\xi_p\Sigma\xi_q\Sigma\xi_F})\Sigma Z}{(\Sigma Z^2-(\Sigma\xi_p+\Sigma\xi_q+\Sigma\xi_F)^2)^2}\\[0.5cm]
    \Sigma\mathcal I_+(\sigma):=\frac{1}{216}\frac{\left(1-\frac{\Sigma m_p\Sigma m_F}{\Sigma\xi_p\Sigma\xi_F}-\frac{\Sigma m_q\Sigma m_F}{\Sigma\xi_q\Sigma\xi_F}+\frac{\Sigma m_p\Sigma m_q}{\Sigma\xi_p\Sigma\xi_q}\right)(\Sigma Z^2+(\Sigma\xi_p+\Sigma\xi_q+\Sigma\xi_F)^2)}{(\Sigma Z^2-(\Sigma\xi_p+\Sigma\xi_q+\Sigma\xi_F)^2)^2}
  \end{array}
\end{equation}
in which
\begin{equation}
  \begin{array}{c}
    \Sigma\xi_p:=\bar\xi\left(\sqrt3+\rho\cos\theta,\omega+\omega\rho\sin\theta\right),\quad
    \Sigma\xi_q:=\bar\xi\left(\sqrt3+\rho\bar r\cos(\theta+\varphi),\omega+\omega\rho\bar r\sin(\theta+\varphi)\right),
    \\[0.3cm]
    \Sigma\xi_F:=\bar\xi\left(\sqrt3+\rho(\cos\theta+\bar r\cos(\theta+\varphi)),\omega+\omega\rho(\sin(\theta)+\bar r\sin(\theta+\varphi))\right),
  \end{array}
\end{equation}
where $\varphi$ is defined in~\-(\ref{eqphi}), and
\begin{equation}
  \begin{array}{>\displaystyle c}
    \bar\xi(\bar k):=\sqrt{\bar m^2(\bar k)+2t_1^2\bar\alpha_1(\bar k)},\quad
    \bar\zeta(\bar k):=2t_2\cos\phi\bar\alpha_1(\bar k),\quad
    \bar m(\bar k):=W-2t_2\sin\phi\bar\alpha_2(\bar k),\\[0.3cm]
    \bar\alpha_1(\bar k_1,\bar k_2):=2\cos\left(\frac\pi3 \bar k_2\right)\left(\cos\left(\frac\pi{\sqrt3}\bar k_1\right)+\cos\left(\frac\pi3 \bar k_2\right)\right)+\frac12,\\[0.5cm]
    \bar\alpha_2(\bar k_1,\bar k_2):=2\sin\left(\frac\pi3 \bar k_2\right)\left(\cos\left(\frac\pi{\sqrt3}\bar k_1\right)-\cos\left(\frac\pi3 \bar k_2\right)\right),
  \end{array}
  \label{eqbarxi}
\end{equation}
\begin{equation}
  \begin{array}{c}
    \Sigma m_p:=\bar m\left(\sqrt3+\rho\cos\theta,\omega+\omega\rho\sin\theta\right),\quad
    \Sigma m_q:=\bar m\left(\sqrt3+\rho\bar r\cos(\theta+\varphi),\omega+\omega\rho\bar r\sin(\theta+\varphi)\right),
    \\[0.3cm]
    \Sigma m_F:=\bar m\left(\sqrt3+\rho(\cos\theta+\bar r\cos(\theta+\varphi)),\omega+\omega\rho(\sin(\theta)+\bar r\sin(\theta+\varphi))\right),
  \end{array}
\end{equation}
and
\begin{equation}
  \Sigma Z:=\Sigma\zeta_p+\Sigma\zeta_q-\Sigma\zeta_F
\end{equation}
with
\begin{equation}
  \begin{array}{c}
    \Sigma\zeta_p:=\bar \zeta\left(\sqrt3+\rho\cos\theta,\omega+\omega\rho\sin\theta\right),\quad
    \Sigma\zeta_q:=\bar \zeta\left(\sqrt3+\rho\bar r\cos(\theta+\varphi),\omega+\omega\rho\bar r\sin(\theta+\varphi)\right),
    \\[0.3cm]
    \Sigma\zeta_F:=\bar \zeta\left(\sqrt3+\rho(\cos\theta+\bar r\cos(\theta+\varphi)),\omega+\omega\rho(\sin(\theta)+\bar r\sin(\theta+\varphi))\right).
  \end{array}
\end{equation}
\bigskip

\indent Let us now check that the functions $\Sigma J\Sigma\mathcal I_\mp$ and $\Sigma f_{\omega,i}$ are smooth. This is almost a direct consequence of lemma~\-\ref{lemmasunrise} and~\-(\ref{eqasymp}), if not for the fact that the sunrise coordinates ignore the periodic nature of the Brillouin zone $\mathcal B$. If $p+q-p_F^\omega$ were equal to $p_F^\omega+(n_1G_++n_2G_-)$ with $G_\pm=(\frac2\pi3,\pm\frac{2\pi}{\sqrt3})$ and $(n_1,n_2)\in\mathbb Z^2\setminus\{(0,0)\}$ then $\mathcal I_\mp(p,q)$ would have a singularity that is not regularized by the sunrise coordinates. However, one readily checks that this cannot happen when $p$ and $q$ are in $\mathcal B_\omega^{(F)}$. All in all, $\Sigma J\Sigma\mathcal I_\mp$ is an analytic function on the closure of the integration domain, and $\Sigma f_{\omega,i}^{(F)}$ is a class-2 Gevrey function.
\bigskip

\makelink{rkphi}{\thepage}
{\bf Remark}: In the discussion above, we assumed that $\mathcal I(p,q)$ is not singular at $p_F^{-\omega}$, which is only true if $\phi\neq0$. If $\phi$ is small, then the derivatives of $\mathcal I(p,q)$ may be very large if one of $p$, $q$ or $p+q-p_F^\omega$ is close to $p_F^{-\omega}$. When $p$ and $q$ are in $\mathcal B_{\omega}^{(F)}$, $p+q-p_F^\omega$ may be arbitrarily close to $p_F^{-\omega}$, which means that $\phi$ must be sufficiently far from $0$ for the accuracy of the computation described above to be good.
\bigskip

\indent Finally, using the $\frac{2\pi}3$ rotation symmetry, we rewrite
\begin{equation}
  A_{F,F}^{(\mp)}=6\int_0^1d\rho\int_{-\frac\pi6}^{\frac\pi2}d\theta\int_{-\frac\pi2}^{\frac\pi2}d\psi\int_0^1dz\ \Sigma f_{\omega,1}^{(F)}(\sigma)\Sigma f_{\omega,2}^{(F)}(\sigma)\Sigma J(\sigma)\Sigma\mathcal I_\mp(\sigma).
\end{equation}
\bigskip

\point The integrand of $A_{R,F}^{(\mp)}$ is only singular if $q=p_F^\omega$ or $p+q-p_F^\omega=p_F^\omega$, because $|p-p_F^\omega|_{\mathcal B}>\frac R2$. We regularize these singularities by switching to polar coordinates corresponding to $q$ and $p+q-p_F^\omega$, which we denote by $(r,\theta,\rho,\varphi)$: if $p+q-p_F^\omega\in\mathcal B_\nu$,
\begin{equation}
  q=p_F^\omega+\omega\frac{2\pi}{3\sqrt3}\rho(\cos\theta,\sin\theta),\quad
  p+q-p_F^\omega=p_F^\nu+\nu\frac{2\pi}{3\sqrt3}r(\cos\varphi,\sin\varphi)
\end{equation}
in terms of which
\begin{equation}
  A_{R,F}^{(\mp)}=\sum_{\nu=\pm}\int_0^{2\pi}d\theta\int_{0}^{2\pi}d\varphi\int_0^{1}dr\int_0^{R(\varphi)}d\rho\ \rho r\Pi f_{\omega}^{(R)}(\varpi)\Pi f_{\omega}^{(F)}(\varpi)\Pi\mathcal I_\mp(\varpi)
\end{equation}
with $\varpi\equiv(r,\theta,\rho,\varphi)$,
\begin{equation}
  \Pi f_{\omega}^{(R)}(\varpi):=\chi_{\frac12}\left(|(\rho\cos\varphi-r\cos\theta,\ \nu-\omega+\nu \rho\sin\varphi-\omega r\sin\theta)|_{\mathbb T}\right),\quad
  \Pi f_{\omega}^{(F)}(\varpi):=\chi_{\frac12}(r)
\end{equation}
where
\begin{equation}
  |\bar k|_{\mathbb T}:=\min\left\{|\bar k+n_1(\sqrt3,1)+n_2(\sqrt3,-1)|,\ (n_1,n_2)\in\mathbb Z^2\right\},
\end{equation}
\begin{equation}
  \begin{array}{>\displaystyle c}
    \Pi\mathcal I_-(\varpi):=\frac{1}{108}\frac{(\Pi\xi_p+\Pi\xi_q+\Pi\xi_F)(\frac{\Pi m_p}{\Pi\xi_p}+\frac{\Pi m_q}{\Pi\xi_q}-\frac{\Pi m_F}{\Pi\xi_F}-\frac{\Pi m_p\Pi m_q\Pi m_F}{\Pi\xi_p\Pi\xi_q\Pi\xi_F})\Pi Z}{(\Pi Z^2-(\Pi\xi_p+\Pi\xi_q+\Pi\xi_F)^2)^2}\\[0.5cm]
    \Pi\mathcal I_+(\varpi):=\frac{1}{216}\frac{\left(1-\frac{\Pi m_p\Pi m_F}{\Pi\xi_p\Pi\xi_F}-\frac{\Pi m_q\Pi m_F}{\Pi\xi_q\Pi\xi_F}+\frac{\Pi m_p\Pi m_q}{\Pi\xi_p\Pi\xi_q}\right)(\Pi Z^2+(\Pi\xi_p+\Pi\xi_q+\Pi\xi_F)^2)}{(\Pi Z^2-(\Pi\xi_p+\Pi\xi_q+\Pi\xi_F)^2)^2}
  \end{array}
\end{equation}
in which
\begin{equation}
  \begin{array}{c}
    \Pi\xi_q:=\bar\xi\left(\sqrt3+r\cos\theta,\omega+\omega r\sin\theta\right),\quad
    \Pi\xi_F:=\bar\xi\left(\sqrt3+\rho\cos\varphi,\nu+\nu \rho\sin\varphi\right),
    \\[0.3cm]
    \Pi\xi_p:=\bar\xi\left(\sqrt3-r\cos\theta+\rho\cos\varphi),\nu-\omega r\sin\theta+\nu \rho\sin\varphi\right),
  \end{array}
\end{equation}
where $\bar\xi$ is defined in~\-(\ref{eqbarxi}),
\begin{equation}
  \begin{array}{c}
    \Pi m_q:=\bar m\left(\sqrt3+r\cos\theta,\omega+\omega r\sin\theta\right),\quad
    \Pi m_F:=\bar m\left(\sqrt3+\rho\cos\varphi,\nu+\nu \rho\sin\varphi\right),
    \\[0.3cm]
    \Pi m_p:=\bar m\left(\sqrt3-r\cos\theta+\rho\cos\varphi),\nu-\omega r\sin\theta+\nu \rho\sin\varphi\right),
  \end{array}
\end{equation}
where $\bar m$ is defined in~\-(\ref{eqbarxi}),
\begin{equation}
  \Pi Z:=\Pi\zeta_p+\Pi\zeta_q-\Pi\zeta_F
\end{equation}
with
\begin{equation}
  \begin{array}{c}
    \Pi\zeta_q:=\bar\zeta\left(\sqrt3+r\cos\theta,\omega+\omega r\sin\theta\right),\quad
    \Pi\zeta_F:=\bar\zeta\left(\sqrt3+\rho\cos\varphi,\nu+\nu \rho\sin\varphi\right),
    \\[0.3cm]
    \Pi\zeta_p:=\bar\zeta\left(\sqrt3-r\cos\theta+\rho\cos\varphi),\nu-\omega r\sin\theta+\nu \rho\sin\varphi\right)
  \end{array}
\end{equation}
where $\bar\zeta$ is defined in~\-(\ref{eqbarxi}), and
\begin{equation}
  R(\theta)=
  \left\{\begin{array}{>\displaystyle l}
    \frac1{\cos(\theta-\frac\pi6)}\quad\mathrm{if\ }\theta\in\left[-\frac\pi6,\frac\pi2\right]\\[0.5cm]
    \frac1{\cos(\theta-\frac{5\pi}6)}\quad\mathrm{if\ }\theta\in\left[\frac\pi2,\frac{7\pi}6\right]\\[0.5cm]
    \frac1{\cos(\theta+\frac\pi2)}\quad\mathrm{if\ }\theta\in\left[\frac{7\pi}6,\frac{11\pi}6\right].
  \end{array}\right.
  \label{eqR}
\end{equation}
Note that $R(\theta)$ is smooth by parts, so, in order to keep the accuracy of the computation high, we must split the integral over $\varphi$:
\begin{equation}
  A_{R,F}^{(\mp)}=3\sum_{\nu=\pm}\int_{0}^{2\pi}d\theta\int_{-\frac\pi6}^{\frac\pi2}d\varphi\int_0^{1}dr\int_0^{R(\varphi)}d\rho\ \rho r\Pi f_{\omega}^{(R)}(\varpi)\Pi f_{\omega}^{(F)}(\varpi)\Pi\mathcal I_\mp(\varpi)
\end{equation}
In which we used the symmetry under $\frac{2\pi}3$ rotations of $p$ and $q$.
\bigskip

\indent By~\-(\ref{eqasymp}) and the fact that $|p-p_F^\omega|_{\mathcal B}>\frac R2$ on the support of $f_\omega^{(R)}$, $\rho r\Pi\mathcal I_\mp$ is an analytic function on the closure of the integration domain, and $\Pi f_{\omega}^{(F)}$ and $\Pi f_{\omega}^{(R)}$ are class-2 Gevrey functions.
\bigskip

\makelink{rkcutoff}{\thepage}
{\bf Remark}: In order to regularize the singularity at $p+q-p_F^\omega=p_F^\omega$, we had to change variables to $(q,p+q-p_F^\omega)$. If, instead of the smooth cutoff function $f_\omega$, we had used a step function, the integration region for $p+q-p_F^\omega$ would have been $\mathcal B$ minus a disk centered around $q$ of radius $R$. This creates trouble, since the parametrization of this disk is singular when $p_F^\omega$ tends to the boundary of the disk. The reason for which we have used a smooth cutoff function is to avoid this problem.
\bigskip

\point The integrand of $A_{R,R}^{(\mp)}$ is only singular if $p+q-p_F^\omega=p_F^\omega$, because $|p-p_F^\omega|_{\mathcal B}>\frac R2$ and $|q-p_F^\omega|_{\mathcal B}>\frac R2$. We regularize this singularity by switching to polar coordinates corresponding to $q$ and $p+q-p_F^\omega$, which we denote by $(r,\theta,\rho,\varphi)$: if $q\in\mathcal B_\eta$ and $p+q-p_F^\omega\in\mathcal B_\nu$, then we define
\begin{equation}
  q=p_F^\eta+\eta\frac{2\pi}{3\sqrt3} r(\cos\theta,\sin\theta),\quad
  p+q-p_F^\omega=p_F^\nu+\nu\frac{2\pi}{3\sqrt3}\rho(\cos\varphi,\sin\varphi)
\end{equation}
in terms of which
\begin{equation}
  A_{R,R}^{(\mp)}=\sum_{\eta,\nu=\pm}\int_0^{2\pi}d\theta\int_{0}^{2\pi}d\varphi\int_{0}^{R(\theta)}dr\int_0^{R(\varphi)}d\rho\ \rho r\Xi f_{\omega,1}^{(R)}(\varpi)\Xi f_{\omega,2}^{(R)}(\varpi)\Xi\mathcal I_\mp(\varpi)
\end{equation}
with $\varpi\equiv(r,\theta,\rho,\varphi)$,
\begin{equation}
  \begin{array}c
    \Xi f_{\omega,1}^{(R)}(\varpi):=\chi_{\frac12}\left(|(\rho\cos\varphi-r\cos\theta,\ \nu-\eta+\nu \rho\sin\varphi-\eta r\sin\theta)|_{\mathbb T}\right),\\[0.3cm]
    \Xi f_{\omega,2}^{(R)}(\varpi):=\chi_{\frac12}\left(|(r\cos\theta,\ \eta-\omega+\eta r\sin\theta)|_{\mathbb T}\right)
  \end{array}
\end{equation}
where
\begin{equation}
  |\bar k|_{\mathbb T}:=\min\left\{|\bar k+n_1(\sqrt3,1)+n_2(\sqrt3,-1)|,\ (n_1,n_2)\in\mathbb Z^2\right\},
\end{equation}
\begin{equation}
  \begin{array}{>\displaystyle c}
    \Xi\mathcal I_-(\varpi):=\frac{1}{108}\frac{(\Xi\xi_p+\Xi\xi_q+\Xi\xi_F)(\frac{\Xi m_p}{\Xi\xi_p}+\frac{\Xi m_q}{\Xi\xi_q}-\frac{\Xi m_F}{\Xi\xi_F}-\frac{\Xi m_p\Xi m_q\Xi m_F}{\Xi\xi_p\Xi\xi_q\Xi\xi_F})\Xi Z}{(\Xi Z^2-(\Xi\xi_p+\Xi\xi_q+\Xi\xi_F)^2)^2}\\[0.5cm]
    \Xi\mathcal I_+(\varpi):=\frac{1}{216}\frac{\left(1-\frac{\Xi m_p\Xi m_F}{\Xi\xi_p\Xi\xi_F}-\frac{\Xi m_q\Xi m_F}{\Xi\xi_q\Xi\xi_F}+\frac{\Xi m_p\Xi m_q}{\Xi\xi_p\Xi\xi_q}\right)(\Xi Z^2+(\Xi\xi_p+\Xi\xi_q+\Xi\xi_F)^2)}{(\Xi Z^2-(\Xi\xi_p+\Xi\xi_q+\Xi\xi_F)^2)^2}
  \end{array}
\end{equation}
in which
\begin{equation}
  \begin{array}{c}
    \Xi\xi_q:=\bar\xi\left(\sqrt3+r\cos\theta,\eta+\eta r\sin\theta\right),\quad
    \Xi\xi_F:=\bar\xi\left(\sqrt3+\rho\cos\varphi,\nu+\nu \rho\sin\varphi\right),
    \\[0.3cm]
    \Xi\xi_p:=\bar\xi\left(\sqrt3-r\cos\theta+\rho\cos\varphi),\nu+\omega-\eta-\eta r\sin\theta+\nu \rho\sin\varphi\right),
  \end{array}
\end{equation}
where $\bar\xi$ is defined in~\-(\ref{eqbarxi}),
\begin{equation}
  \begin{array}{c}
    \Xi m_q:=\bar m\left(\sqrt3+r\cos\theta,\eta+\eta r\sin\theta\right),\quad
    \Xi m_F:=\bar m\left(\sqrt3+\rho\cos\varphi,\nu+\nu \rho\sin\varphi\right),
    \\[0.3cm]
    \Xi m_p:=\bar m\left(\sqrt3-r\cos\theta+\rho\cos\varphi),\nu+\omega-\eta-\eta r\sin\theta+\nu \rho\sin\varphi\right),
  \end{array}
\end{equation}
where $\bar m$ is defined in~\-(\ref{eqbarxi}),
\begin{equation}
  \Xi Z:=\Xi\zeta_p+\Xi\zeta_q-\Xi\zeta_F
\end{equation}
with
\begin{equation}
  \begin{array}{c}
    \Xi\zeta_q:=\bar\zeta\left(\sqrt3+r\cos\theta,\eta+\eta r\sin\theta\right),\quad
    \Xi\zeta_F:=\bar\zeta\left(\sqrt3+\rho\cos\varphi,\nu+\nu \rho\sin\varphi\right),
    \\[0.3cm]
    \Xi\zeta_p:=\bar\zeta\left(\sqrt3-r\cos\theta+\rho\cos\varphi),\nu+\omega-\eta-\eta r\sin\theta+\nu \rho\sin\varphi\right)
  \end{array}
\end{equation}
where $\bar\zeta$ is defined in~\-(\ref{eqbarxi}), and $R$ is defined in~\-(\ref{eqR}). Here, again, since $R(\theta)$ is only smooth by parts, we must split the integral over $\theta$ and $\varphi$:
\begin{equation}
  A_{R,R}^{(\mp)}=3\sum_{\eta,\nu=\pm}\sum_{a=0,1,2}\int_{(4a-1)\frac\pi6}^{(4a+3)\frac\pi6}d\theta\int_{-\frac\pi6}^{\frac\pi2}d\varphi\int_{0}^{R(\theta)}dr\int_0^{R(\varphi)}d\rho\ \rho r\Xi f_{\omega,1}^{(R)}(\varpi)\Xi f_{\omega,2}^{(R)}(\varpi)\Xi\mathcal I_\mp(\varpi)
\end{equation}
In which we used the symmetry under $\frac{2\pi}3$ rotations of $p$ and $q$.
\bigskip

\indent By~\-(\ref{eqasymp}) and the fact that $|p-p_F^\omega|_{\mathcal B}>\frac R2$ and $|q-p_F^\omega|_{\mathcal B}>\frac R2$ on the support of $f_\omega^{(R)}$, $r\Xi\mathcal I_\mp$ is an analytic function on the closure of the integration domain, and $\Xi f_{\omega,1}^{(R)}$ and $\Xi f_{\omega,2}^{(R)}$ are class-2 Gevrey functions.
\bigskip

\subsection{Strategy of the numerical computation}
\indent The integrations are carried out using Gauss-Legendre quadratures (see section~\-\ref{secintegrationgl}).

\subsection{Usage and examples}
\indent We will now describe some basic usage cases of {\tt hhtop z1-z2}. For a full description of the options of {\tt hhtop}, see the {\tt man} page.
\subseqskip

\subsubsection{Basic usage}
\indent The value of the parameters can be set via the {\tt -p} flag. Here is an example\par
\medskip
\indent{\tt hhtop z1-z2 -p "omega=1;t1=1.;t2=.1;phi=1;"}\par
\indent{\tt hhtop z1+z2 -p "omega=1;t1=1.;t2=.1;phi=1;"}\par
\medskip
The parameters that are not specified by the {\tt-p} flag are set to their default value: $\omega=1$, $t_1=1$, $t_2=0.1$, $\phi=\frac\pi2$, $W=\omega3\sqrt{3}t_2\sin\phi$.

\subsubsection{Precision of the computation}
\indent The precision of the computation can be controlled by three parameters: the precision of the numbers manipulated by {\tt hhtop} (set via the {\tt -P} flag, see section~\-\ref{secprecision}), the order of the integration, and the tolerance of the computation of abcissa and weights.
\bigskip

\point{\bf Order of the integration.} The order of the integration, that is, the value of the number $N$ introduced in section~\-\ref{secintegrationgl}, can be specified via the {\tt -O} flag. Its default value is 10. The difference of the value of the integral at different orders is a good measure of the numerical error. Example:\par
\medskip
\indent{\tt hhtop z1-z2 -O 30}\par
\indent{\tt hhtop z1+z2 -O 30}
\bigskip

\point{\bf Tolerance of the abcissa and weights.} A Newton scheme is used to compute the abcissa and weights of the Gauss-Legendre integration. The scheme halts when the difference $|x_{n+1}-x_n|$ (see section~\-\ref{secnewton}) is smaller than a number, called the {\it tolerance} of the algorithm, or when the toal number of steps exceeds a given threshold. The tolerance can be set via the {\tt-t} flag, and the maximal number of steps via the {\tt -N} flag. Their default values are $10^{-11}$ and 1000000. If the tolerance or the maximal number of steps are too small, and the precision of multi-precision floats is too low, then the iteration may not converge. Example:\par
\medskip
\indent{\tt hhtop z1-z2 -t 1e-30 -N 2000000 -O 100 -P 256}\par
\indent{\tt hhtop z1+z2 -t 1e-30 -N 2000000 -O 100 -P 256}

\subsubsection{Using double precision floats instead of multi-precision floats}
\indent Using the {\tt -D} command-line flag, {\tt hhtop} can be instructed to use {\tt long double}'s instead of MPFR floats. Whereas one then loses the ability of adjusting the precision, the computation time can be drastically reduced. Example:\par
\medskip
\indent{\tt hhtop -D z1-z2 -p "sinphi=1.;"}\par
\indent{\tt hhtop -D z1+z2 -p "sinphi=1.;"}\par
\bigskip
The precision of {\tt long double}s is compiler dependent, see section~\-\ref{secprecision}.

\vfill\eject


\section{Algorithms}
\indent In this section, we describe the algorithms used by {\tt hhtop}. Their implementation is provided by the {\tt libinum} library.
\subseqskip

\subsection{Newton scheme}\label{secnewton}
\indent The Newton algorithm is used to compute roots: given a real function $f$ and an initial guess $x_0$ for the root, the Newton scheme produces a sequence $(x_n)$:
\begin{equation}\begin{array}{>\displaystyle c}
x_{n+1}:=x_n-\frac{f(x_n)}{\partial_x f(x_n)}
\end{array}\label{eqnewton}\end{equation}
which, provided the sequence converges, it tends to a solution of $f(x)=0$, with a quadratic rate of convergence
\begin{equation}
|x_{n+1}-x_{n}|\leqslant c_n|x_{n}-x_{n-1}|^2
\label{eqconvnewton}\end{equation}
where
\begin{equation}
c_n:=\frac12\frac{\mAthop{\displaystyle\mathrm{sup}}_{x\in[x_{n+1},x_{n}]}|\partial^2_x f(x)|}{\mAthop{\displaystyle\mathrm{inf}}_{x\in[x_{n+1},x_n]}|\partial_xf(x)|}.
\label{eqboundNewton}\end{equation}

\subsection{Gauss-Legendre integration}\label{secintegrationgl}
\indent The Gauss-Legendre method allows us to compute
\begin{equation}
\int_{-1}^1dx\ f(x)
\label{eqgenericint}\end{equation}
for $f:[-1,1]\to\mathbb{R}$. Having fixed an {\it order} $N\in\mathbb N\setminus\{0\}$, let $\{x_1,\cdots,x_N\}$ denote the set of roots of the $N$-th Legendre polynomial $P_N$,and let
\begin{equation}
w_i=\frac2{(1-x_i^2)P_N'(x_i)}
\label{eqweights}\end{equation}
for $i\in\{1,\cdots,N\}$. One can show that, if $f$ is a polynomial of order $\leqslant2N-1$, then
\begin{equation}
\int_{-1}^1dx\ f(x)=\sum_{i=1}^Nw_if(x_i).
\label{eqgenericint}\end{equation}
\bigskip

\indent If $f$ is an analytic function, then one can show that the error decays exponentially as $N\to\infty$. However, in the computation of $z_1-z_2$ (see section~\-\ref{secz1z2}), we use Gauss-Legendre quadratures to integrate a class-2 Gevrey function, so we will need to generalize this result. Let us first define class-$s$ Gevrey functions on $[-1,1]$, as $\mathcal C^\infty$ functions, for which there exist $C_0,C>0$, such that $\forall n\in\mathbb N$,
\begin{equation}
  \mathop{\mathrm{sup}}_{x\in[-1,1]}\left|\frac{d^nf(x)}{dx^n}\right|\leqslant C_0C^n(n!)^s.
\end{equation}
Note that the set of analytic functions on $[-1,1]$ is equal to the set of class-1 Gevrey functions on $[-1,1]$. We assume that $s\in\mathbb N\setminus\{0\}$.
\bigskip

\indent The basic strategy to estimate the error
\begin{equation}
  E_N(f):=\left|\int_{-1}^1dx\ f(x)-\sum_{i=1}^Nw_if(x)\right|
\end{equation}
is to approximate $f$ using Chebyshev polynomials, bound the error of this approximation for Gevrey functions, and use an estimate of the error when $f$ is a Chebyshev polynomial. This is done in detail in appendix~\-\ref{appGL}, and we find~\-(see lemma~\-\ref{lemmaGL})
\begin{equation}
  E_N(f)\leqslant c_0c_1^{s-1}(2N)^{1-\frac1s}e^{-b(2N)^{\frac1s}}s!.
\end{equation}


\subsection{Precision}\label{secprecision}
\indent The numerical values manipulated by {\tt hhtop} are represented as multi-precision floats (using the GNU MPFR library). The number of bits allocated to each number, that is, the number of digits used in the computation, can be specified using the {\tt -P} flag. The default precision is 53 bits. Example:\par
\medskip
\indent{\tt hhtop phase -P 128}\par
\bigskip

\indent This behavior can be changed using the {\tt -D} flag, in which case the numerical values are represented as {\tt long double}, which have a fixed precision, but yield faster computation times. Example:\par
\medskip
\indent{\tt hhtop -D phase}\par
\bigskip
The precision of {\tt long double}'s is compiler-dependent, and can be checked using the {\tt -Vv} flag:
\medskip
\indent{\tt hhtop -Vv}\par
\bigskip
Using the GNU GCC compiler, version 5.3.0, on the {\tt x86-64} architecture, the precision of {\tt long double}'s is 64.


\vfill\eject

\appendix
\section{Proof of (\expandonce{\ref{eqineqab}})}\label{appintk0}
\indent In this appendix, we show that~\-(\ref{eqcondt}) holds, then~\-(\ref{eqineqab}) does as well. To alleviate the notation, we will drop the $_{t_1,t_2,W,\phi}$ indices as well as the `$(k)$'. We have
\begin{equation}
\xi^2-\zeta^2=m^2+t_1^2|\Omega|^2-4t_2^2\cos^2\phi\alpha_1^2
\label{eqximinuszeta1}\end{equation}
which, using $|\Omega|^2=2\alpha_1$, becomes
\begin{equation}
\xi^2-\zeta^2=m^2+2\alpha_1(t_1^2-2t_2^2\alpha_1).
\label{eqximinuszeta2}\end{equation}
Furthermore, $0\leqslant\alpha_1\leqslant\frac92$ (both $0$ and $\frac92$ are reached, respectively at 0 and $p_F:=(\frac{2\pi}3,\frac{2\pi}{3\sqrt3})$). This implies that $\xi^2>\zeta^2$.

\section{Sunrise coordinates}\label{appsunrise}
\indent In this appendix, we discuss the {\it sunrise} coordinates, which are used to compute sunrise Feynman diagrams. Such diagrams give rise to an integral of the form
\begin{equation}
  \int dpdq\ F(p,q,|p|,|q|,|p+q|)
\end{equation}
where $\rho r F(p,q,\rho,r,\gamma)$ is an analytic function of $p$, $q$, $\rho$, $r$ and $\gamma$, and
\begin{equation}
  |(p_1,p_2)|:=\sqrt{p_1^2+p_2^2}.
\end{equation}
However, since $|p|$, $|q|$ and $|p+q|$ are not analytic, the derivatives of $F(p,q,|p|,|q|,|p+q|)$ are, typically, unbounded, which can cause the error in the numerical evaluation of the integral uncontrollably large. In order to avoid this problem, we introduce coordinates, $(\rho,\theta,\psi,z)$, called {\it sunrise coordinated}, which are such that $p$, $q$, $|p|$, $|q|$, $|p+q|$, as well as the Jacobian of the change of variables, are analytic functions of $(\rho,\theta,\psi,z)$. Expressed using the sunrise coordinates, the integral of $F$ can be computed with good numerical accuracy.
\bigskip

{\bf Remark}: Note that if, instead of the sunrise coordinates, one used the (simpler) polar coordinates $p=\rho(\cos\theta,\sin\theta)$ and $q=r(\cos\varphi,\sin\varphi)$, then $|p+q|=\sqrt{\rho^2+r^2+2\rho r\cos(\theta-\varphi)}$, which has a divergent second derivative at $(\rho,\theta)=(r,-\varphi)$. Polar coordinates, therefore, do not do the trick.
\bigskip

{\bf Remark}: The sunrise coordinates are introduced in the following lemma, which is only stated for the case $|p|>|q|$. The integration over the regime $|q|<|p|$ can be performed by exchanging $p$ and $q$.
\bigskip

\theo{Lemma}\label{lemmasunrise}
  Let $\mathcal B_R:=\{p\in\mathbb R^2,\ |p|<R\}$. We define the map $\mathcal S$
  \begin{equation}
    \begin{array}{rrcl}
      \mathcal S:&
      \{(p,q)\in\mathcal B_R^2,\ |p|>|q|\}
      &\longrightarrow&
      (0,R)\times[0,2\pi)\times[-\frac\pi2,\frac\pi2]\times(0,1)\\[0.3cm]
      &(p,q)&\longmapsto&(\rho,\theta,\psi,z)
    \end{array}
  \end{equation}
  with
  \begin{equation}
    \rho:=|p|\in(0,R),
    \label{eqrho}
  \end{equation}
  $\theta\in[0,2\pi)$ is the unique solution of
  \begin{equation}
    p=\rho(\cos\theta,\sin\theta),
    \label{eqtheta}
  \end{equation}
  if $\varphi$ denotes the angle between $p$ and $q$, then $\psi\in[-\frac\pi2,\frac\pi2]$ is the unique solution of
  \begin{equation}
    \cos\psi=\sqrt{\frac{|p+q|-|p|+|q|}{2|q|}},\quad
    \mathrm{sign}(\psi)=\mathrm{sign}(\sin\varphi),
    \label{eqpsi}
  \end{equation}
  and
  \begin{equation}
    z:=1-\frac{1-\sqrt{1-\frac{|q|}\rho\sin^2\psi}}{1-\cos\psi}\in(0,1).
    \label{eqz}
  \end{equation}
  The map $\mathcal S$ is invertible, its inverse is analytic, and is such that, if $(p,q)=\mathcal S^{-1}(\rho,\theta,\psi,z)$, then $|p|$, $|q|$ and $|p+q|$ are analytic functions of $(\rho,\theta,\psi,z)$. Furthermore, the Jacobian
  \begin{equation}
    J:=\left|\det\left(\frac{\partial(p_1,p_2,q_1,q_2)}{\partial(\rho,\theta,\psi,z)}\right)\right|
  \end{equation}
  is an analytic function of $(\rho,\theta,\psi,z)$. In addition, $\mathcal S^{-1}$, $|p|$, $|q|$, $|p+q|$ and $J$, as functions of $(\rho,\theta,\psi,z)$, can be continued analytically to $[0,R]\times[0,2\pi)\times[-\frac\pi2,\frac\pi2]\times[0,1]$. Explicitly,
  \begin{equation}
    p_1=\rho\cos\theta,\quad
    p_2=\rho\sin\theta,\quad
    q_1=\rho\bar r\cos(\theta+\varphi),\quad
    q_2=\rho\bar r\sin(\theta+\varphi),
    \label{eqpq}
  \end{equation}
  \begin{equation}
    |p|=\rho,\quad
    |q|=\rho\bar r,\quad
    |p+q|=\rho(1+\bar r\cos(2\psi))
    \label{eqnorms}
  \end{equation}
  and
  \begin{equation}
    J=4\rho^3\frac{\bar r(1+\bar r\cos(2\psi))}{(1+\cos\psi)\sqrt{1+\bar r\cos^2\psi}}
    \label{eqjacobian}
  \end{equation}
  with
  \begin{equation}
    \bar r:=(1-z)(1+zh(\psi)),\quad
    h(\psi):=\frac{1-\cos\psi}{1+\cos\psi},\quad
    t:=1-(1-z)(1-\cos\psi),
    \label{eqr}
  \end{equation}
  and
  \begin{equation}
    \cos\varphi:=\cos(2\psi)-\frac{\bar r}2\sin^2(2\psi),\quad
    \sin\varphi:=t\sin(2\psi)\sqrt{1+\bar r\cos^2\psi}.
    \label{eqphi}
  \end{equation}
\endtheo
\bigskip

\indent\underline{Proof}: In order to prove the lemma, we will compose several changes of coordinates. The {\it sunrise} coordinates described above are obtained by combining these intermediate changes of variables.
\bigskip

\point The first, consists in changing $p$ to polar coordinates, which yields~\-(\ref{eqrho}), (\ref{eqtheta}) and the first two equations of~\-(\ref{eqpq}), and contributes a factor $\rho$ to the Jacobian:
\begin{equation}
  \int_{\mathcal B_R}dp\int_{\mathcal B_{|p|}}dq\ F=\int_0^R d\rho\int_0^{2\pi} d\theta\int_{\mathcal B_\rho} dq\ \rho\ F.
\end{equation}
\bigskip

\point We then change variables to
\begin{equation}
  (\rho,\theta,q_1,q_2)\longmapsto(\rho,\theta,r,\gamma)
\end{equation}
with
\begin{equation}
  r:=|q|,\quad
  \gamma:=|p+q|=\sqrt{\rho^2+r^2+2\rho(q_1\cos\theta+q_2\sin\theta)}
\end{equation}
so that
\begin{equation}
  \int_{\mathcal B_R}dp\int_{\mathcal B_{|p|}}dq\ F= \int_0^R d\rho \int_0^{2\pi}d\theta \int_0^\rho dr \int_{\rho-r}^{\rho+r}d\gamma\ \frac \gamma{|\sin\varphi|}\ F
\end{equation}
where $\varphi$ is the angle between $p$ and $q$:
\begin{equation}
  \cos\varphi=\frac{\gamma^2-\rho^2-r^2}{2r\rho},\quad
  |\sin\varphi|=\sqrt{1-\cos^2\varphi}=\frac1{2r\rho}\sqrt{4r^2\rho^2-(\gamma^2-r^2-\rho^2)^2}
\end{equation}
which we rewrite as
\begin{equation}
  |\sin\varphi|=\frac1{2r\rho}
  \sqrt{(\rho+r+\gamma)(\rho-r+\gamma)(-\rho+r+\gamma)(\rho+r-\gamma)}.
\end{equation}
\bigskip

\point We then adimensionalize $r$ and $\gamma$, that is, we change to $\bar r,\bar \gamma$ in such a way that $\bar r,\bar \gamma\in(0,1)$:
\begin{equation}
  \bar r:=\frac r\rho,\quad
  \bar \gamma:=\frac{\gamma-\rho+r}{2r}
\end{equation}
in terms of which
\begin{equation}
  \int_{\mathcal B_R}dp\int_{\mathcal B_{|p|}}dq\ F= \int_0^R d\rho \int_0^{2\pi}d\theta \int_0^1 d\bar r \int_0^1d\bar \gamma\ \frac{2\rho^3\bar r(1-\bar r+2\bar r\bar \gamma)}{|\sin\varphi|}\ F
\end{equation}
and
\begin{equation}
  |\sin\varphi|=2\sqrt{\bar \gamma(1-\bar \gamma)(1-\bar r+\bar r\bar \gamma)(1+\bar r\bar \gamma)}.
  \label{eqphi2}
\end{equation}
\bigskip

\point At this point, the singularities have all been shifted to $|\sin\varphi|$: $p$, $q$, $|p|$, $|q|$ and $|p+q|$ are analytic functions of $\rho$, $\bar r$, $\bar \gamma$, $\cos\theta$, $\sin\theta$, $\cos\varphi$ and $\sin\varphi$, and the only one of these that is singular is $\sin\varphi$, because of the square root in~\-(\ref{eqphi2}). We first note that $\sqrt{1+\bar r\bar \gamma}>1$, so that factor is not singular. In order to regularize the divergence in the other terms, we change variables to
\begin{equation}
  \cos\psi:=\sqrt{\bar \gamma},\quad
  \sin\psi:=\mathrm{sign}(\sin\varphi)\sqrt{1-\bar \gamma},\quad
  t:=\sqrt{1-\bar r(1-\bar \gamma)}
\end{equation}
after which
\begin{equation}
  \int_{\mathcal B_R}dp\int_{\mathcal B_{|p|}}dq\ F= \int_0^R d\rho \int_0^{2\pi}d\theta \int_{-\frac\pi2}^{\frac\pi2} d\psi \int_{\cos\psi}^1dt\ 4\rho^3\frac{(1-t^2)}{\sin^4\psi}\frac{\left(1+\frac{(1-t^2)}{\sin^2\psi}(2\cos^2\psi-1)\right)}{\sqrt{1+\frac{(1-t^2)}{\sin^2\psi}\cos^2\psi}}\ F.
\end{equation}
\bigskip

\point Finally, we adimensionalize $t$:
\begin{equation}
  z:=1-\frac{1-t}{1-\cos\psi}
\end{equation}
so that
\begin{equation}
  \begin{largearray}
    \int_{\mathcal B_R}dp\int_{\mathcal B_{|p|}}dq\ F= \int_0^R d\rho \int_0^{2\pi}d\theta \int_{-\frac\pi2}^{\frac\pi2} d\psi \int_0^1dz
    \\[0.3cm]\hfill
    4\rho^3\frac{(1-z)(1+zh(\psi))}{1+\cos\psi}\frac{(1+(1-z)(1+zh(\psi))(2\cos^2\psi-1))}{\sqrt{1+(1-z)(1+zh(\psi))\cos^2\psi}}\ F.
  \end{largearray}
  \label{eqintfinal}
\end{equation}
\bigskip

\point Equations~\-(\ref{eqpq}) through~\-(\ref{eqphi}) follow from~\-(\ref{eqintfinal}). The analyticity of $p$, $q$, $|p|$, $|q|$, $|p+q|$ and $J$ is a simple comsequence of~\-(\ref{eqpq}), (\ref{eqnorms}) and~\-(\ref{eqjacobian}).\qed

\vfill\eject


\section{Estimate of the error of Gauss-Legendre quadratures for Gevrey functions}\label{appGL}
\indent In this appendix, we compute the error of Gauss-Legendre quadratures when used to integrate class-$s$ Gevrey functions. A class-$s$ Gevrey function on $[-1,1]$ is a $\mathcal C^\infty$ function that satisfies, $\forall n\in\mathbb N$,
\begin{equation}
  \mathop{\mathrm{sup}}_{x\in[-1,1]}\left|\frac{d^nf(x)}{dx^n}\right|\leqslant C_0C^n(n!)^s.
\end{equation}

\bigskip

\theo{Lemma}\label{lemmaGL}
  Let $f$ be a class-$s$ Gevrey function with $s\in\mathbb N\setminus\{0\}$. There exist $c_0,c_1,b>0$, which are independent of $s$, and $N_0>0$, which is independent of $s$ and $f$, such that, if $N\geqslant N_0$, then
  \begin{equation}
    E_N(f)\leqslant c_0c_1^{s-1}(2N)^{1-\frac1s}e^{-b(2N)^{\frac1s}}s!.
  \end{equation}
  In particular, if $f$ is analytic (i.e. $s=1$), then
  \begin{equation}
    E_N(f)\leqslant c_0e^{-2bN}.
  \end{equation}
\endtheo
\bigskip

\indent\underline{Proof}:\par\penalty10000
\medskip\penalty10000
\point We approximate $f$ by Chebyshev polynomials:
\begin{equation}
  f(x)=\frac{c_0}2+\sum_{j=1}^{\infty}c_jT_j(x)
  \label{eqcheby}
\end{equation}
where $T_j$ is the $j$-th Chebyshev polynomial:
\begin{equation}
  T_j(x):=\cos(j\arccos(x)),\quad
  c_j:=\frac2\pi\int_0^\pi d\theta\ f(\cos\theta)\cos(j\theta).
\end{equation}
Note that~\-(\ref{eqcheby}) is nothing other than the Fourier cosine series expansion of $F(\theta):=f(\cos(\theta))$, which is an even, periodic, class-$s$ Gevrey function on $[-\pi,\pi]$, whose $j$-th Fourier coefficient for $j\in\mathbb Z$ is equal to $\frac12c_{|j|}$. Furthermore, using a well-known estimate of the decay of Fourier coefficients of class-$s$ Gevrey functions (see e.g.~\-[\cite{Ta87}, Theorem~\-3.3]), there exists $b_0,b>0$ such that
\begin{equation}
  c_j\leqslant b_0e^{-bj^{\frac1s}}.
  \label{eqcjbound}
\end{equation}
\bigskip

\point Furthermore, since order-$N$ Gauss-Legendre quadratures are exact on polynomials of order $\leqslant 2N-1$, we have, formally,
\begin{equation}
  E_N(f)=\sum_{j=2N}^\infty c_jE_N(T_j).
\end{equation}
As was proved by A.R.~\-Curtis and P.~Rabinowitz~\-[\cite{CR72}], if $N$ is large enough, then
\begin{equation}
  E_N(T_j)\leqslant\pi
\end{equation}
which, by~\-(\ref{eqcjbound}), implies that
\begin{equation}
  E_N(f)\leqslant\pi\sum_{j=2N}^\infty c_j\leqslant\pi b_0\sum_{j=2N}^\infty e^{-bj^{\frac1s}}.
\end{equation}
Furthermore, if $\nu_{N,s}^s:=\lfloor(2N)^{\frac1s}\rfloor^s$ denotes the largest integer that is $\leqslant 2N$ and has an integer $s$-th root, then
\begin{equation}
  \sum_{j=2N}^\infty e^{-bj^{\frac1s}}\leqslant
  \sum_{j=\nu_{N,s}^s}^\infty e^{-bj^{\frac1s}}\leqslant
  \sum_{k=\nu_{N,s}}^\infty(k^s-(k-1)^s)e^{-bk}\leqslant
  s\sum_{k=\nu_{N,s}}^\infty k^{s-1}e^{-bk}.
\end{equation}
We then estimate
\begin{equation}
  \sum_{k=\nu_{N,s}}^\infty k^{s-1}e^{-bk}=\frac{d^{s-1}}{d(-b)^{s-1}}\sum_{k=\nu_{N,s}}^\infty e^{-bk}
  \leqslant (s-1)!\left(\nu_{N,s}+\frac1{1-e^{-b}}\right)^{s-1}\frac{e^{-b\nu_{N,s}}}{1-e^{-b}}
\end{equation}
which concludes the proof of the lemma.\qed

\vfill\eject

\references
\BBlography

\end{document}