\documentclass[sort&compress]{elsarticle}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{color}
\usepackage{natbib}
\usepackage{hyperref}

\newcommand{\dd}{\mathrm{d}}
\newcommand{\im}{\mathrm{Im}}

\newcommand{\LMQ}{l_q}
\newcommand{\ann}[1]{{\tiny\color{red}PM: #1}}
\newcommand{\annAM}[1]{{\tiny\color{blue}AM: #1}}
\begin{document}
% \begin{flushleft}
% TTP13-0XX\\
% SFB/CPP-13-XX\\
% arXiv:YYMM.NNNN
% \end{flushleft}
\begin{frontmatter}
\title{The QED vacuum polarization function at four loops and the
  anomalous magnetic moment at five loops} 
\author[Moskau]{P.~Baikov}
\author[Muenchen]{A.~Maier}
\author[Zeuthen]{P.~Marquard} \address[Moskau]{Skobeltsyn Institute of Nuclear Physics, Lomonosov Moscow State University,
1(2), Leninskie gory, Moscow 119991, Russian Federation}
\address[Muenchen]{Physik
Department T31, James-Franck-Stra\ss{}e, Technische Universit\"at M\"unchen, D–85748 Garching, Germany} 
% \address[Karlsruhe]{Institut f\"ur Theoretische
%   Teilchenphysik, Karlsruhe~Institute~of~Technology~(KIT), 76128
%   Karlsruhe, Germany}
\address[Zeuthen]{Deutsches Elektronen Synchrotron DESY, Platanenallee 6, D15637 Zeuthen, Germany}
% \begin{flushleft}
% TTP13-0XX\\
% SFB/CPP-13-XX\\
% arXiv:YYMM.NNNN
% \end{flushleft}
\begin{abstract}
  The anomalous moment of the muon is one of the most fundamental
  observables. It has been measured experimentally with a very high
  precision and on theory side the contributions from perturbative QED
  have been calculated up to five-loop level by numerical
  methods. Contributions to the muon anomalous magnetic moment from
  certain diagram classes are also accessible by alternative methods. In
  this paper we present the evaluation of contributions to the QCD
  corrections due to insertions of the vacuum polarization function at
  five-loop level.
\end{abstract}
\begin{keyword}
Perturbative calculations, QED, lepton anomalous moment 
%\PACS
\end{keyword}
\end{frontmatter}

\section{Introduction}
\label{sec:intro}
The anomalous magnetic moments of electron and muon are among the best
experimentally measured quantities and are also very well understood from
theoretical side. In particular, the perturbative QED corrections have
been under consideration for a very long time. The leading-order QED
corrections have been considered by Schwinger in
Ref.~\cite{Schwinger:1948iu}, the next-to-leading order corrections not
much later in Refs.~\cite{Petermann:1957hs,Sommerfield:1957zz}. At three
loops, i.e. next-to-next-to-leading order, the QED corrections have been
first calculated numerically in Ref.~\cite{Kinoshita:1995ym} and later
analytically in
Refs.~\cite{Laporta:1996mq,Melnikov:2000qh,Marquard:2007uj}.  At four and
five loops the complete QED contributions have been calculated
numerically in Refs.~\cite{Aoyama:2012wk,Aoyama:2012wj} while only
partial analytical results
exist~\cite{Lautrup:1974ic,Kinoshita:1990ur,Kawai:1991wd,Faustov:1990zs,Broadhurst:1992za,Baikov:1995ui,Laporta:1993ds,Baikov:2012rr,Lee:2013sx}. For
a more thorough review of the current status see
e.g. Ref.~\cite{Jegerlehner:2009ry}.

In the case of the muon anomalous moment, next to the universal
contribution, contributions from diagrams with closed electron loops are
of particular interest. A subgroup of this class of diagrams,
contributions from corrections to the photon propagator, can be easily
obtained using already available building blocks
\cite{Lautrup:1974ic,Barbieri:1974nc}. At five loops the leading
contributions have been obtained from considering the asymptotic form~\cite{Baikov:2012rr}, which resulted in unexpected discrepancies with
the results obtained in Ref.\cite{Aoyama:2010zp}. The aim of this paper
is to improve the predictions for the anomalous magnetic moment of the
muon made in Ref.~\cite{Baikov:2012rr} and to resolve the
discrepancies. To this extend we approximately reconstruct the photon
vacuum polarization at four loops using all available information.

This paper is organized as follows: In
Section~\ref{sec:qed-vacu-polar} we collect the necessary information
to construct an approximation for the vacuum polarization function at
four loops. In Section~\ref{sec:muon} we present the calculation of the
contribution to the anomalous magnetic muon from vacuum polarization
insertions. 

\section{QED vacuum polarization at four loops}
\label{sec:qed-vacu-polar}
We define the vacuum polarization function $\Pi(q^2)$ as usual by
\begin{equation}
  \label{eq:4}
  (q^\mu q^\nu - q^2 g^{\mu\nu}) \Pi(q^2) = i \int \dd x \,\langle 0 | e^{i q x} T j^\mu(x) j^\nu(0) | 0 \rangle   \,,
\end{equation}
with the current $j^\mu = \bar\psi \gamma^\mu \psi$ and write it as an
expansion in the fine-structure constant $\alpha$
\begin{equation}
  \label{eq:5}
  \Pi(q^2) = \frac{\alpha}{\pi}\Pi^{(1)}(q^2) 
+\left(\frac{\alpha}{\pi}\right )^2\Pi^{(2)}(q^2) 
+\left(\frac{\alpha}{\pi}\right )^3\Pi^{(3)}(q^2) 
+\left(\frac{\alpha}{\pi}\right )^4\Pi^{(4)}(q^2) 
+{\cal O}(\alpha^5)  \,.
\end{equation}
At three and four loops the vacuum polarization can only be calculated
in a low-energy, high-energy or threshold expansion at the moment.  The
available information can then be used to construct an approximate
function using Pad\'e approximation. In the following we will present
the results for the corresponding expansions and construct an
approximating function. The renormalization of lepton mass and
fine-structure constant $\alpha$ is performed in the on-shell
scheme. Note that in this scheme the condition $\Pi(0) = 0$ is being
imposed.

\subsection{Low-energy expansion}
\label{sec:low-energy-expansion}
In the low-energy limit the polarization function can be expanded in a
power series in $z=q^2 / ( 4m_q^2 ) < 1$
\begin{equation}
  \label{eq:6}
  \Pi^{(n)}_{\mathrm{le}} = \frac{1}{16\pi^2}\sum_{k=1} \Pi^{(n)}_{\mathrm{le},k} z^k \,.
\end{equation}
For the QED case at hand, the result for the non-singlet
contribution can in principle be obtained from the QCD results given in
Ref.~\cite{Maier:2009fz}. But since the results are not given expressed
through colour factors they cannot easily be translated to the case of
QED. The low-energy expansion was therefore recalculated for the case at
hand. The calculation follows a well-established path. The Feynman
diagrams are generated using {\tt qgraf}~\cite{Nogueira:1993} and mapped
onto six topologies using {\tt q2e} and {\tt
exp}~\cite{Seidensticker:1999bb,Harlander:1997zb}. Then a {\tt
FORM}~\cite{vermaseren-form} program is used to apply projectors and
take traces. The resulting scalar integrals are then reduced to master
integrals using {\tt Crusher}~\cite{crusher}, which implements Laporta's
algorithm~\cite{Laporta:2001dd} for solving integration-by-parts
identities. The needed master integrals have been calculated in
Refs.~\cite{Chetyrkin:2006dh,Schroder:2005hy,Laporta:2002pg,Chetyrkin:2004fq,Kniehl:2005yc,Schroder:2005db,Bejdakic:2006vg,Kniehl:2006bf,Kniehl:2006bg}. Combining
all steps and performing the renormalization of all quantities in the
on-shell scheme, the first three moments are obtained { \small
  \begin{equation}
    \label{eq:234}
    \begin{split}
&\Pi^{(4)}_{\mathrm{le},1} =
\Bigg(-\frac{30853 a_1^4}{22680}+\frac{30853 \pi ^2 
a_1^2}{22680}-\frac{64 \pi ^2 a_1}{135}-\frac{30853 
a_4}{945}-\frac{77255063 \zeta_3}{4762800}  \\
&+\frac{473237 \pi 
^4}{1360800}+\frac{866 \pi ^2}{1215}-\frac{3619201}{529200}\Bigg) 
n_l^2+\Bigg(\frac{3376 a_1^5}{2835}-\frac{2178299 
a_1^4}{544320}-\frac{3376 \pi ^2 a_1^3}{1701} \\
&+\frac{2178299 \pi ^2 
a_1^2}{544320}-\frac{12191 \pi ^4 a_1}{8505}-\frac{1922 \pi ^2 
a_1}{405}-\frac{2178299 a_4}{22680}-\frac{27008 a_5}{189}+\frac{\pi ^2 \zeta
_3}{30}\\
&-\frac{69694097 \zeta_3}{1209600}+\frac{33931 \zeta_5}{189}+\frac{66975707 \pi ^4}{65318400}+\frac{7297 \pi 
^2}{3240}+\frac{3562169}{2721600}\Bigg) n_l+ \\
&\mathbf{si} \Bigg(-\frac{739 
a_1^4}{1080}+\frac{739 \pi ^2 a_1^2}{1080}-\frac{739 
a_4}{45}-\frac{2017831 \zeta_3}{213840}+\frac{35 \pi 
^4}{216}+\frac{664837}{641520}\Bigg) n_l^2 \\
&+\Bigg(-\frac{1291 \zeta
_3}{1944}-\frac{32 \pi ^2}{2025}+\frac{83971}{78732}\Bigg) n_l^3 \,,
\end{split}
\end{equation}
% 
%
%
\begin{equation}
\begin{split}
&\Pi^{(4)}_{\mathrm{le},2} =
\Bigg(-\frac{130829911 a_1^4}{8709120}+\frac{130829911 \pi ^2 
a_1^2}{8709120}-\frac{128 \pi ^2 a_1}{315}-\frac{130829911 
a_4}{362880}\\
&-\frac{5763324918049 \zeta_3}{26824089600}+\frac{2208846791 \pi 
^4}{522547200}+\frac{10846 \pi 
^2}{14175}-\frac{42466787908001}{1086375628800}\Bigg) 
n_l^2 \\
&+\Bigg(\frac{19593724 a_1^5}{2027025}-\frac{60530131639 
a_1^4}{24908083200}-\frac{19593724 \pi ^2 
a_1^3}{1216215}+\frac{60530131639 \pi ^2 
a_1^2}{24908083200}\\
&-\frac{72540947 \pi ^4 a_1}{6081075}-\frac{3058 \pi ^2 
a_1}{675}-\frac{60530131639 a_4}{1037836800}-\frac{156749792 
a_5}{135135}+\frac{\pi ^2 \zeta_3}{35}\\
&-\frac{112244692092317 \zeta
_3}{6974263296000}+\frac{14315837 \zeta_5}{10010}-\frac{427149037853 \pi 
^4}{2988969984000}+\frac{839 \pi 
^2}{378}+\frac{94330906317547}{3487131648000}\Bigg) n_l\\
&+\mathbf{si} 
\Bigg(\frac{97011619 a_1^4}{174182400}-\frac{97011619 \pi ^2 
a_1^2}{174182400}+\frac{97011619 a_4}{7257600}+\frac{796232393699 \zeta
_3}{92990177280}-\frac{745372259 \pi 
^4}{4180377600} \\
&+\frac{5881974201847}{2092278988800}\Bigg) 
n_l^2+\Bigg(-\frac{19669747 \zeta_3}{21288960}-\frac{64 \pi 
^2}{4725}+\frac{3284183491}{2586608640}\Bigg) n_l^3 \,,
\end{split}
\end{equation}
% 
%
\begin{equation}
\begin{split}
&\Pi^{(4)}_{\mathrm{le},3} =
\Bigg(-\frac{1875259367 a_1^4}{17740800}+\frac{1875259367 \pi ^2 
a_1^2}{17740800}-\frac{1024 \pi ^2 a_1}{2835}-\frac{1875259367 
a_4}{739200}\\
&-\frac{377287031234107 \zeta_3}{245887488000}+\frac{95566793477 
\pi ^4}{3193344000}+\frac{1175348 \pi 
^2}{1488375}-\frac{2388270016962373}{9958443264000}\Bigg) 
n_l^2 \\
&+\Bigg(\frac{12135758968 a_1^5}{172297125}-\frac{245658223193 
a_1^4}{6616209600}-\frac{12135758968 \pi ^2 
a_1^3}{103378275}+\frac{245658223193 \pi ^2 
a_1^2}{6616209600}\\
&-\frac{45718624634 \pi ^4 a_1}{516891375}-\frac{720316 
\pi ^2 a_1}{165375}-\frac{245658223193 a_4}{275675400}-\frac{97086071744 
a_5}{11486475}+\frac{8 \pi ^2 \zeta_3}{315}\\
&-\frac{1531450738927589 \zeta
_3}{3368252160000}+\frac{119140260224 \zeta
_5}{11486475}+\frac{176076905389817 \pi ^4}{31757806080000}+\frac{5337 \pi 
^2}{2450}\\
&+\frac{300002162759308069}{1500556337280000}\Bigg) n_l+\mathbf{si} 
\Bigg(\frac{22845879073 a_1^4}{6096384000}-\frac{22845879073 \pi ^2 
a_1^2}{6096384000}\\
&+\frac{22845879073 a_4}{254016000}+\frac{343009147408727 
\zeta_3}{6246309888000}-\frac{787819133821 \pi 
^4}{731566080000}+\frac{7118016595194017}{758926651392000}\Bigg) 
n_l^2\\
&+\Bigg(-\frac{7731286469 \zeta_3}{6227020800}-\frac{512 \pi 
^2}{42525}+\frac{81866930683}{50438868480}\Bigg) n_l^3 \,,
    \end{split}
  \end{equation}
}
where $n_l$ denotes contributions from closed lepton loops and where
we used the abbreviations $a_1 = \log 2$ and $ a_n =
\mathrm{Li}_n(1/2)$. $\mathbf{si}$ marks contributions from singlet diagrams.
\subsection{High-energy expansion}
\label{sec:high-energy-expans}
In the high-energy region we write the result in the form
\begin{equation}
  \label{eq:7}
  \Pi^{(n)}_{\mathrm{he}} = \frac{1}{16\pi^2}\sum_{k=0} \Pi^{(n)}_{\mathrm{he},k} z^{-k}  \,.
\end{equation}
The coefficients of Eq.~(\ref{eq:7}) can be expressed through
four-loop massless propagator integrals which, in turn,
can be reduced to 28 master integrals.
This reduction has been done by evaluating sufficiently many
terms of the $1/D$ expansion
\cite{Baikov:2005nv}
of the corresponding coefficient functions
\cite{Baikov:1996rk}.
The master integrals are known analytically from
\cite{Baikov:2010hf}.

As the result, the leading two coefficients of Eq.~(\ref{eq:7}) are
{\small
\begin{equation}
\begin{split}
& \Pi^{(4)}_{\mathrm{he},0}  =\mathbf{si}  \Bigg(\frac{73 a_1^4}{144}-\frac{73}{144} \pi ^2 
a_1^2+\frac{73 a_4}{6}-\frac{2 \zeta _3^2}{3}+\frac{5309 \zeta 
_3}{1120}+\frac{5 \zeta _5}{3}+\frac{2}{3} \zeta _3 \log (-4 
z) \\
&-\frac{11}{36} \log (-4 z)-\frac{2237 \pi 
^4}{17280}+\frac{1963}{3780}\Bigg)n_l^2 +n_l^2 \Bigg(\frac{53 
a_1^4}{60}-\frac{53}{60} \pi ^2 a_1^2+\frac{16 \pi ^2 a_1}{27}+\frac{106 
a_4}{5}-\zeta _3^2 \\
&+\frac{29129 \zeta _3}{1800}-\frac{125 \zeta 
_5}{18}-\frac{19}{12} \zeta _3 \log (-4 z)+\frac{5}{3} \zeta _5 \log (-4 
z)-\frac{1}{48} \log ^2(-4 z)-\frac{1}{12} \log (-4 z)\\
&-\frac{2161 \pi 
^4}{10800}-\frac{179 \pi ^2}{324}+\frac{3361}{900}\Bigg)+n_l 
\Bigg(-\frac{32 a_1^5}{225}+\frac{1559 a_1^4}{1080}+\frac{32}{135} \pi 
^2 a_1^3-\frac{1559 \pi ^2 a_1^2}{1080}\\
&+\frac{106 \pi ^4 
a_1}{675}+\frac{59 \pi ^2 a_1}{12}+\frac{1559 a_4}{45}+\frac{256 
a_5}{15}-\frac{\pi ^2 \zeta _3}{24}+\frac{6559 \zeta _3}{320}-\frac{1603 
\zeta _5}{120}-\frac{35 \zeta _7}{4}\\
&+\frac{23}{128} \log (-4 z)-\frac{59801 
\pi ^4}{129600}-\frac{157 \pi ^2}{72}-\frac{71189}{34560}\Bigg)+n_l^3 
\Bigg(-\frac{15109 \zeta _3}{22680}-\frac{5 \zeta _5}{9}\\
&-\frac{1}{9} \zeta 
_3 \log ^2(-4 z)+\frac{19}{27} \zeta _3 \log (-4 z)-\frac{1}{108} \log ^3(-4 
z)+\frac{11}{72} \log ^2(-4 z)-\frac{151}{162} \log (-4 z)\\
&+\frac{8 \pi 
^2}{405}+\frac{75259}{68040}\Bigg) \,,
\end{split}
\end{equation}
\begin{equation}
\begin{split}
& \Pi^{(4)}_{\mathrm{he},1}  =n_l^2 \Bigg(\frac{a_1^4}{9}-\frac{1}{9} \pi ^2 a_1^2+\frac{8 \pi ^2 
a_1}{9}+\frac{8 a_4}{3}-\frac{\zeta _3^2}{6}+\frac{3023 \zeta 
_3}{432}-\frac{415 \zeta _5}{54}\\
&-\frac{167}{72} \zeta _3 \log (-4 
z)-\frac{35}{36} \zeta _5 \log (-4 z)+\frac{3}{16} \log ^3(-4 
z)-\frac{15}{16} \log ^2(-4 z)+\frac{1}{4} \pi ^2 \log (-4 
z)\\
&+\frac{539}{144} \log (-4 z)+\frac{5 \pi ^4}{108}-\frac{103 \pi 
^2}{108}-\frac{491}{216}\Bigg)+n_l \Bigg(-\frac{3}{4} \pi ^2 a_1 \log (-4 
z)-\frac{a_1^4}{2}+\frac{1}{2} \pi ^2 a_1^2\\
&+\frac{31 \pi ^2 a_1}{4}-12 
a_4-2 \zeta _3^2-\frac{\pi ^2 \zeta _3}{16}-\frac{3647 \zeta 
_3}{192}-\frac{1825 \zeta _5}{192}+\frac{10073 \zeta _7}{384}+\frac{41}{32} 
\zeta _3 \log (-4 z)\\
&+\frac{35}{16} \zeta _5 \log (-4 z)-\frac{9}{32} \log 
^3(-4 z)+\frac{9}{64} \log ^2(-4 z)+\frac{15}{32} \pi ^2 \log (-4 
z)-\frac{685}{128} \log (-4 z)\\
&-\frac{\pi ^4}{48}-\frac{673 \pi 
^2}{192}-\frac{47}{36}\Bigg)+\mathbf{si}\Bigg(-\frac{9 \zeta _3^2}{2}-\frac{13 \zeta 
_3}{3}-\frac{5 \zeta _5}{12}+\frac{147 \zeta _7}{16}+\frac{1}{3}\Bigg) 
 n_l^2\\
&+n_l^3 \Bigg(-\frac{4 \zeta _3}{9}+\frac{2}{3} \zeta _3 
\log (-4 z)-\frac{1}{36} \log ^3(-4 z)+\frac{13}{72} \log ^2(-4 
z)-\frac{317}{216} \log (-4 z)\\
&+\frac{4 \pi ^2}{135}+\frac{25}{27}\Bigg) \,,
\end{split}
\end{equation}
}
where $n_l$ denotes contributions from closed lepton loops and where
we used the abbreviations $a_1 = \log 2$ and $ a_n =
\mathrm{Li}_n(1/2)$.  $\mathbf{si}$ marks contributions from singlet diagrams.

\subsection{Threshold expansion}
\label{sec:threshold-expansion}
The polarization function in the threshold region can be written as
\begin{equation}
  \label{eq:pi_thr}
  \Pi^{(n)}_{\mathrm{thr}} = \sum_{k=2-n} \Pi^{(n)}_{\mathrm{thr},k} \sqrt{1-z\,}^{\,k}  \,.
\end{equation}
In non-relativistic quantum field theory the NNLO expression for the 
threshold cross section is known for an arbitrary $SU(N)$ gauge
group~\cite{Hoang:2000yr}. This means that the derivation of the threshold expansion
is essentially the same as for the QCD case, which is described in 
Refs.~\cite{Hoang:2008qy,Kiyo:2009gb}. The only additional steps consist
of setting the group invariants to their $U(1)$ values and converting
the coupling to the on-shell scheme.

The resulting coefficients read
\begin{align}
  \label{pi3_thr_coeff}
  \Pi^{(3)}_{\text{thr},-1} =& \frac{\pi^3}{384}\*n_l\,,\\
  \Pi^{(3)}_{\text{thr},0} =& \bigg(\frac{17}{768}+\frac{11}{128}\*a_1\bigg)\*n_l\*\log(1-z)+\frac{11}{512}\*n_l\*\log^2(1-z)+\text{const}\,,\\
\Pi^{(3)}_{\text{thr},1} =
&\bigg(\frac{17}{192\*\pi}-\frac{397}{2304}\*\pi+\frac{5}{24}\*\pi\*a_1-\frac{13}{2304}\*\pi^3+\frac{41}{64\*\pi}\*\zeta_3\bigg)\*n_l\notag\\
&+\frac{11}{192}\*\pi\*n_l\*\log(1-z)
+\bigg(-\frac{11}{72\*\pi}+\frac{\pi}{72}\bigg)\*n_l^2\,,
\end{align}
for $n=3$ and
\begin{align}
  \label{pi4_thr_coeff_-2}
  \Pi^{(4)}_{\text{thr},-2} =& \frac{1}{128}\*\pi^2\*\zeta_3\*n_l\,,\\
  \label{pi4_thr_coeff_-1}
  \Pi^{(4)}_{\text{thr},-1} =&
  \bigg(\frac{7}{2304}\*\pi^3-\frac{11}{384}\*\pi^3\*a_1-\frac{11}{64}\*\pi\*\zeta_3\bigg)\*n_l-\frac{11}{768}\*\pi^3\*n_l\*\log(1-z)\,,\\
  \label{pi4_thr_coeff_0}
  \Pi^{(4)}_{\text{thr},0} =&
\bigg(-\frac{467}{18432}+\frac{115}{384}\*a_1-\frac{121}{256}\*a_1^2-\frac{2207}{9216}\*\pi^2+\frac{5}{48}\*\pi^2\*a_1\notag\\
&+\frac{9}{4096}\*\pi^4+\frac{131\*\zeta_3}{512}\bigg)\*n_l\*\log(1-z)\notag\\
&+\bigg(-\frac{839}{9216}+\frac{\pi^2}{144}\bigg)\*n_l^2\*\log(1-z)\notag\\
&+\bigg(\frac{115}{1536}-\frac{121}{512}\*a_1+\frac{11}{768}\*\pi^2\bigg)\*n_l\*\log^2(1-z)\notag\\
&-\frac{121}{3072}\*n_l\*\log^3(1-z)+\text{const}\,,
\end{align}
for $n=4$, where we again used $a_1 = \log 2$.

\subsection{Pad\'e approximation}
\label{sec:pade-approximation}
Having all building blocks at hand the polarization function can be
reconstructed using Pad\'e approximation~\cite{Baker:1975,Broadhurst:1993mw,Fleischer:1994ef,Broadhurst:1994qj,Baikov:1995ui,Chetyrkin:1995ii,Hoang:2008qy,Kiyo:2009gb}. For
the four-loop contribution we closely follow the procedure outlined in
Ref.~\cite{Kiyo:2009gb}. In the three-loop case we introduce slight
modifications to accommodate the large amount of information in the
available low- and high-energy expansions. We first give a brief review
of the method as used in the four-loop case and then discuss the changes
for the application to the three-loop contribution.

\subsubsection{Approximation procedure at four loops}
\label{sec:pade-approximation-4}
As in Ref.~\cite{Kiyo:2009gb}, we first split $\Pi^{(4)}(z)$ into two parts,
\begin{equation}
  \label{eq:pi_split}
  \Pi^{(4)}(z) = \Pi^{(4)}_{\text{reg}}(z) + \Pi^{(4)}_{\text{log}}(z)\,,
\end{equation}
using the ansatz
\begin{align}
\label{eq:pi4_log}
  \Pi^{(4)}_{\text{log}}(z) =& \sum_{j=-1}^0\sum_{i=1}^{3+2\*j} k_{ij} {\cal K}^{(4)}_{ij}(z) +
  \sum_{n=0}^{1}\sum_{m=0}^{3} d_{mn} {\cal D}^{(4)}_{mn}(z)\,,\\
\label{eq:pi4_log_thr}
  {\cal K}^{(4)}_{ij}(z) =& \Pi^{(2)}(z)^i\*G(z)^j \times
  \begin{cases}
    \left(A_1+\frac{1}{z}\right)&\text{for }i=3,j=0\\
    \left(A_2+\frac{1}{z}\right)& \text{for }i=1,j=-1\\
    \quad1&\text{otherwise}
  \end{cases}\,,
\\
\label{eq:pi4_log_he}
  {\cal D}^{(4)}_{mn}(z) =&
  \big(z\,\*G(z)\big)^m\*\left(1-1/z\right)^{\lceil\frac{m}{2}\rceil}\frac{1}{z^n}\times
  \begin{cases}
    \left(1+\frac{1}{B_1z}\right)& \text{for }n=3\\
    \left(1+\frac{1}{B_2z}\right)& \text{for }n=2\\
    \quad1&\text{otherwise}
  \end{cases}\,,
\end{align}
with 
%arbitrary auxiliary parameters $A_1, A_2, B_1, B_2$, 
the known two-loop polarization function $\Pi^{(2)}(z)$ as
introduced in Eq.~(\ref{eq:5}) and
\begin{equation}
  \label{eq:G_def}
  G(z) = \frac{1}{2z}\frac{\log\left(u\right)}{\sqrt{1-\frac{1}{z}}}\,,
\qquad u = \frac{\sqrt{1-\frac{1}{z}}-1}{\sqrt{1-\frac{1}{z}}+1}\,.
\end{equation}

The coefficients $k_{ij}$ can be fixed so that all known
logarithms in the threshold expansion
(Eqs.~\eqref{eq:pi_thr},~\eqref{pi4_thr_coeff_-2}--~\eqref{pi4_thr_coeff_0})
 are absorbed into the first term on the right-hand side of
 Eq.~\eqref{eq:pi4_log}. Expanding Eq.~\eqref{eq:pi4_log_thr} in the
 threshold region generates the required logarithms and half-integer
 powers of $1-z$, as can be seen from the expansions
\begin{align}
  \label{eq:pi_G_thr}
  G(z) &= \frac{\pi}{2}\frac{1}{\sqrt{1-z}}+{\cal O}((1-z)^0)\,,\\
\Pi^{(2)}(z)=&-\frac{1}{16}\log\left(1-z\right)+\text{const} +{\cal O}(\sqrt{1-z}) \,.
\end{align}

In a similar way the factor $z\,\*G(z)$ in Eq.~\eqref{eq:pi4_log_he}
generates logarithms in the high-energy region:
\begin{equation}
  \label{G_he}
  z\,\*G(z) = -\frac{1}{2}\log(-4\*z)+{\cal O}\left(\frac{1}{z}\right)\,.
\end{equation}
This means we can choose the coefficients $d_{mn}$ with $m>0$ in such a way that
these logarithms are also absorbed into $\Pi^{(4)}_{\text{log}}(z)$. The
remaining coefficients $d_{0n}$ are fixed by requiring that
$\Pi^{(4)}_{\text{reg}}(z)$ has no poles for $z\to 0$ and, more specifically,
$\Pi^{(4)}_{\text{reg}}(0) = 0$.

To estimate the error of the
approximation we vary the parameters $A_i,B_i$ independently with
\begin{align}
  \label{eq:pade_params}
  A_i \in & \{-1 \pm 1,\, -1\pm 4,\, -1 \pm 16,\, -1 \pm 64\}\,,\notag\\
  B_i \in & \{\pm1,\, \pm 4,\, \pm 16,\, \pm64\}\,.
\end{align}

In a second step we define
\begin{equation}
  \label{eq:pade_aux_4l}
  P_1(\omega) =
  \left(\frac{1-\omega}{1+\omega}\right)^2\*\left(\Pi^{(4)}_{\text{reg}}\big(z(w)\big)
    - \Pi^{(4)}_{\text{reg}}(-\infty)\right)
\end{equation}
with $z(\omega) = 4\*\omega/(1+\omega)^2$ and construct Pad\'e
approximants
\begin{equation}
  \label{eq:pade_def}
  p_{n,m}(\omega) =  \frac{\sum_{i=0}^{n-1} a_i \omega^i+\omega^n}{\sum_{i=0}^m b_i \omega^i}
\end{equation}
using the constraints
\begin{align}
  \label{eq:pade_cons}
  p_{n,m}(1)=&P_1(1), & p'_{n,m}(1)=&P'_1(1),\notag\\
  p_{n,m}(-1)=&P_1(-1), & p^{(i)}_{n,m}(0)=&P^{(i)}_1(0),& &i=0,\dots,3 \,,
\end{align}
and requiring the absence of terms proportional to $z^{-1/2}$ in the
high-energy expansion of $\Pi^{(4)}_{\text{reg}}$. In total there are
eight constraints, i.e. we obtain Pad\'e approximants with $n+m=7$.

If $p_{n,m}(\omega)$ has poles inside the unit circle $|\omega|<1$ the
corresponding reconstructed polarization function shows unphysical
singularities in the complex $z$ plane. We therefore discard such
approximants. Furthermore, we require $|p_{n,m}(\omega)| < 1.5$ for
$|\omega|=1$ in order to remove approximants with pronounced additional
peaks above the physical threshold at $z=1$. A notable effect of this
cut is the elimination of all Taylor approximants $p_{7,0}(\omega)$.

Since the approximants are only available in
numerical form we show the general features in Fig.~\ref{fig:app}.
\begin{figure}
  \centering
  \includegraphics[width=0.8\textwidth]{plotPi3.eps}
  \caption{Pad\'e approximation for $\Pi^{(4)}\left (-\frac{x^2}{1-x}m_\mu^2\right )$. We show in the top
    half the approximants and in the bottom half the relative error
    with respect to the local mean of all approximants
    obtained.}
  \label{fig:app}
\end{figure}
At the top of the figure we show the envelope of all Pad\'e
approximants that have been calculated and compare with the low- and
high-energy expansions. In the bottom half we show the maximal deviation
from the mean in percent. As can be seen there are only significant
deviations in the range $x>0.8$, which correspond to the uncertainty
due to a limited number of terms in the high-energy expansion.

\subsubsection{Modifications at three loops}
\label{sec:pade-approximation-3}
At three loops much deeper expansions can be used for the construction of
the Pad\'e approximants and solving the resulting large system of
constraints becomes more involved. Our strategy will be to incorporate as
much of the threshold information as possible into the logarithmic
function $\Pi^{(3)}_{\text{log}}(z)$, so that all constraints are
imposed at $\omega=0$. The resulting system can be solved very
efficiently using well-established techniques of one-point Pad\'e
approximation.

As in the four-loop case we first split the polarization function into
two parts, using
\begin{align}
  \label{eq:pi3_log}
  \Pi^{(3)}_{\text{log}}(z) =& \sum_{i\geq0,j} k_{ij} {\cal K}^{(3)}_{ij}(z) +
  \sum_{m,n} d_{mn} {\cal D}^{(3)}_{mn}(z) \,,\\
  {\cal K}^{(3)}_{ij}(z) =& \Pi^{(2)}(z)^i\*G(z)^j \times
  \begin{cases}
    (1-z)^j & \text{for }j>0\\
    \quad 1 & \text{otherwise}
  \end{cases}\,,
\\
  {\cal D}^{(3)}_{ij}(z) =&
  \big(z\,\*G(z)\big)^m\*\left(1-1/z\right)^{\lceil\frac{m}{2}\rceil}\frac{1}{z^n}
\end{align}
to absorb the logarithms and the threshold singularity. Since we expect
to obtain sufficiently many different Pad\'e approximants for a reliable error
estimate we refrain from introducing additional parameters at this point.

To map all available information from the low- and high-energy expansions onto
$\omega=0$, we define
\begin{align}
  \label{eq:pade_aux_3l}
  P_{30}(\omega) &= z(\omega)^{31}\*\left(\Pi^{(3)}_{\text{reg}}\big(z(\omega)\big)
    - \sum_{i=0}^{30}\frac{H^{(3)}_i}{z(\omega)^i}\right)\,,\notag\\
 H^{(3)}_i &= \frac{1}{i!}\left(\frac{\partial}{\partial(1/z)}\right)^i\Pi^{(3)}_{\text{reg}}(z)\Bigg|_{z\to -\infty}\,.
\end{align}
The Taylor approximant is fixed by imposing
\begin{equation}
  \label{eq:pade_taylor}
  p^{(i)}_{61,0}(0) = P^{(i)}_{30}(0)\,,\qquad i=0,\dots,61\,.
\end{equation}
We can deduce all further approximants with $n+m=61$ and $n+m=60$
using Baker's recursion formula~\cite{baker:1970}\footnote{Note that
the denominator in the first equation of Eq.~(13) 
in Ref.~\cite{baker:1970} contains an obvious typo. The proper expression is
$\bar{\eta}_{2j-1}\*\theta_{2j-2}(x) - x\*\bar{\eta}_{2j-2}\*\theta_{2j-1}(x)$.
}
\begin{align}
  \label{eq:pade_rec}
p_{n-\frac{j}{2},\frac{j}{2}} 
&= \frac{\hat{\eta}_j}{\hat\theta_j}
= \frac{\hat\eta_{j-2} -
  \omega\*\,\hat\eta_{j-1}}{\hat\theta_{j-2} -
  \omega\*\,\hat\theta_{j-1}} & j \text{ even}\,,\\
p_{n-\left\lceil\frac{j}{2}\right\rceil,\left\lfloor\frac{j}{2}\right\rfloor} 
&= \frac{\hat{\eta}_j}{\hat\theta_j}
= \frac{\hat\eta_{j-2} -
  \hat\eta_{j-1}}{\hat\theta_{j-2} -
  \hat\theta_{j-1}} & j \text{ odd}\,,
\end{align}
where the numerator $\hat{\eta}_j$ of the Pad\'e
approximant is normalised as in Eq.~\eqref{eq:pade_def}.

It should be noted that the known threshold constant at order $1-z$ is
not included in the approximants. We find however that the information
from the threshold region has virtually no effect on the reconstructed
polarization functions in the Euclidean regime $z<0$. This is contrary
to the four-loop case, where the approximants are far less constrained
by the low- and high-energy expansions.

\section{The anomalous magnetic moment of the muon at five loops}
\label{sec:muon}
The QED corrections to the anomalous magnetic moment can be calculated
in perturbation theory and can thus be written in form of a power series
in the fine structure constant $\alpha$
\begin{equation}
  a_\mu = \sum_{k=1}^\infty \left ( \frac{\alpha}{\pi}\right )^k a_\mu
  ^{(2 k)} \,,
\end{equation}
where $a_\mu^{(2 k)}$ can be further decomposed -- following the
conventions in Ref.~\cite{Aoyama:2012wk} -- as 
\begin{equation}
  a_\mu^{(2 k)} = A_1^{(2 k)}  
  + A_2^{(2 k)}(m_e / m_\mu) + A_2^{(2 k)}(m_\tau / m_\mu)  
  + A_3^{(2 k)}(m_e / m_\mu,m_\tau / m_\mu) \,.
\end{equation}
$A_1^{(2 k)}$ contains the universal contributions, which in case of
the muon anomalous magnetic moment only contain muon loops. The
diagrams contributing to $A_2^{(2 k)}(m_e / m_\mu)$ and $A_2^{(2
  k)}(m_\tau / m_\mu)$ have at least one electron or tau loop,
respectively. In $A_3^{(2 k)}(m_e / m_\mu,m_\tau / m_\mu)$
contributions from diagrams with both electron and tau loops are
collected. In this paper we are mainly interested in contributions to 
$A_2^{(2 k)}(m_e / m_\mu)$ without any muon loops. 

The contributions to the anomalous magnetic moment of the
muon due to photon polarization effects can be calculated
(cp. Fig.~\ref{fig:1}) by using
\begin{equation}
\label{eq:1}
a_\mu = \frac{\alpha}{\pi} \int_{0}^1 \dd x (1-x)  \left [ - \Pi
  (s_x) \right ] \,, \, s_x = - \frac{x^2}{1-x} m_\mu^2 \,.
\end{equation}
Multiple insertions of one-particle irreducible vacuum polarization can
be obtained by using
\begin{equation}
\label{eq:2}
a_\mu = \frac{\alpha}{\pi} \int_{0}^1 \dd x (1-x)  \left [ - \Pi^{\mathrm{1PI}}
  (s_x) \right ]^n \,,
\end{equation}
leading to e.g. at two-loop order
\begin{equation}
\label{eq:3}
a_\mu = \frac{\alpha}{\pi} \int_{0}^1 \dd x (1-x)  \left [ -
  \Pi^{\mathrm{2L}}
  (s_x) + 2 \Pi^{\mathrm{1L}}
  (s_x) \Pi^{\mathrm{1L}}
  (s_x)\right ] \,.
\end{equation}
\begin{figure}
  \centering
  \includegraphics[width=0.3\linewidth]{proto}
  \caption{Prototype diagram}
  \label{fig:1}
\end{figure}
The classes of diagrams accessible by this method are shown in
Fig.~\ref{fig:2}. 
\begin{figure}
  \centering
  I(a)\includegraphics[width=0.15\linewidth]{classes}
  I(b)\includegraphics[width=0.15\linewidth]{classes2}
  I(c)\includegraphics[width=0.15\linewidth]{classes5}
  I(d)\includegraphics[width=0.15\linewidth]{classes6}
  I(e)\includegraphics[width=0.15\linewidth]{classes3}
  I(f)\includegraphics[width=0.15\linewidth]{classes7}
  I(g)\includegraphics[width=0.15\linewidth]{classes8}
  I(h)\includegraphics[width=0.15\linewidth]{classes9}
  I(i)\includegraphics[width=0.15\linewidth]{classes4}
  I(j)\includegraphics[width=0.15\linewidth]{classes10}
  \caption{Classes of diagrams accessible by the used method.}
  \label{fig:2}
\end{figure}
\newcommand{\zwf}[3]{\ensuremath{ #1 {{\scriptstyle + #2} \atop
      {\scriptstyle - #3} } }}

\begin{table}
  \centering
\begin{tabular}{l||l|l|ll||}
  &this work&Ref.~\cite{Baikov:2012rr}&
  \multicolumn{2}{|l||}{Refs.~\cite{Kinoshita:2005sm,Aoyama:2008hz,Aoyama:2010zp,Aoyama:2008gy}}\\\hline\hline
  I(a)&20.142 813&20.183 2&20.142 93(23)&\cite{Kinoshita:2005sm}\\\hline
  I(b)&27.690 061&27.718 8&27.690 38(30)&\cite{Kinoshita:2005sm}\\\hline
  I(c) &\phantom{0}4.742 149&\phantom{0}4.817 59&\phantom{0}4.742 12(14)&\cite{Kinoshita:2005sm}\\\hline
  I(d)+I(e)&\phantom{0}6.241 470&\phantom{0}6.117 77&\phantom{0}6.243 32(101)(70)&\cite{Kinoshita:2005sm}\\\hline
  I(e)&-1.211 249&-1.331 41&-1.208 41(70)&\cite{Kinoshita:2005sm}\\\hline
  I(f)+I(g)+I(h)&\zwf{\mbox{\phantom{0}4.446 9}}{6}{4}&\phantom{0}4.391 31&\phantom{0}4.446 68(9)(23)(59)&\cite{Kinoshita:2005sm,Aoyama:2008hz}\\\hline
  I(i)&\zwf{\mbox{\phantom{0}0.074 6}}{8}{19}&\phantom{0}0.252 37&\phantom{0}0.087 1(59)&\cite{Aoyama:2010zp}\\\hline
  I(j)&\zwf{\mbox{-1.246 9}}{4}{3}&-1.214 29&-1.247 26(12)&\cite{Aoyama:2008gy}
\end{tabular}
%   \begin{tabular}{l||l|l|l||}
%     &this work&Ref.~\cite{Baikov:2012rr}&Refs. cited in ~\cite{Baikov:2012rr}\\\hline\hline
% I(a)&20.142 813&20.183 2&20.142 93(23)\\\hline
% I(b)&27.690 061&27.718 8&27.690 38(30)\\\hline
% I(c) &4.742 149&4.817 59&4.742 12(14)\\\hline
% I(d)+I(e)&6.241 470&6.117 77&6.243 32\\\hline
% I(e)&-1.211 249&-1.331 41&-1.208 41(70)\\\hline
% I(f)+I(g)+I(h)&\zwf{\mbox{4.446 9}}{6}{4}&4.391 31&4.446 68\\\hline
% I(i)&\zwf{\mbox{0.074 6}}{8}{19}&0.252 37&0.087 1(59)\\\hline
% I(j)&\zwf{\mbox{-1.246 9}}{4}{3}&-1.214 29&-1.247 26(12)
%   \end{tabular}
  \caption{Results for $A_2^{(10)}(m_e/m_\mu)$ with pure electronic
    insertions. Please note that the authors of Ref.~\cite{Baikov:2012rr} only used the asymptotic form of $\Pi(s)$ and did not provide any error estimate.}
  \label{tab:Table1}
\end{table}
\begin{table}
  \centering
  \begin{tabular}{c||l|l||}
    &this work&Ref.~\cite{Aoyama:2012wk}\\\hline
I(a)&22.566 976&22.566 973 (3)
\\\hline
I(b)&30.667 093&30.667 091 (3)\\\hline
I(c)&\phantom{0}5.141 395&\phantom{0}5.141 395 (1)\\\hline
I(e)&-0.931 839& -0.931 2 (24)
  \end{tabular}
  \caption{Results for $A_2^{(10)}(m_e/m_\mu)$ including electronic and
    muonic contributions.}
  \label{tab:Table2}
\end{table}
\begin{table}
  \centering
  \begin{tabular}{l||l|l||}
    &this work&Ref.~\cite{Aoyama:2012wj}\\\hline\hline
I(a)&0.000 471&0.000 470 94 (6)
\\\hline
I(b)&0.007 010&0.007 010 8 (7)
\\\hline
I(c) &0.023 467&0.023 468 (2)
\\\hline
I(d)+I(e)&0.014 093&0.014 098(5)(4)
\\\hline
I(e)&0.010 291&0.010 296 (4)
\\\hline
I(f)+I(g)+I(h)&0.037 795&0.037 833(20)(6)(13)
\\\hline
I(i)&0.017 167&0.017 47 (11)
\\\hline
I(j)&0.000 390&0.000 397 5 (18)
\end{tabular}
  \caption{Results for the universal contributions $A_1^{(10)}$.}
  \label{tab:Table3}
\end{table}
In the following, we compare the results obtained in our analysis with
results from previous calculations. The numbers shown are obtained by
numerically integrating over the best available approximation. In case
there are several equivalent approximations the result is obtained by
taking the mean of all values obtained. The errors are then calculated
by taking the difference between the mean and the smallest and largest
values obtained, respectively.

For classes I(a)--I(j) we present our results for the case of purely
electronic vacuum polarization insertions in Tab.~\ref{tab:Table1}. We
compare our results with the values obtained in
Ref.~\cite{Baikov:2012rr}, which relies only on the leading
asymptotics of $\Pi(q^2)$, and results obtained by purely numerical
calculations.

Since the results for classes I(a)-I(c) are obtained by numerical
integrating the exact analytical expression for the one- and two-loop
vacuum polarization they are exact. 

In the case of classes I(d) and I(e) the used three-loop Pad\'es are
highly constrained by a large number of terms in the low- and
high-energy expansion. Thus the error from the spread between different
approximants is negligible. As can be seen the results are in good
agreement with the results obtained in Ref.~\cite{Kinoshita:2005sm} and
the analysis in Ref.~\cite{Baikov:2012rr} can clearly be improved by
including sub-leading contributions.

In case of classes I(f)-I(j) the used four-loop Pad\'es are less
precise but also here we find good agreement within the quoted errors
with the results obtained in Ref.~\cite{Aoyama:2012wk}. In all cases
one finds a significant improvement when comparing to the leading
logarithmic approximation used in Ref.~\cite{Baikov:2012rr}.

For classes I(a)-(c) and I(e) we can obtain the full result for
$A_2^{(10)}(m_e/m_\mu)$ including muonic contributions. These results are
presented in Tab.~\ref{tab:Table2}. In Tab.~\ref{tab:Table3} we
present our results for the universal corrections and compare with the
results given in Ref.~\cite{Aoyama:2012wj}. In both cases the
discussion as for the purely electronic contributions can essentially
be repeated and also here good agreement with results available in the
literature is observed.

As a check of our setup we repeated the analysis of
Ref.~\cite{Baikov:1995ui} and found complete agreement with their work.

\section{Conclusions}
\label{sec:conclusions}
We presented results for a certain set of five-loop contributions to
the anomalous magnetic moment of the muon that stem from corrections
to the vacuum polarization of the photon. We have shown that an
improved treatment of the vacuum polarization, including more than its
asymptotic form, leads to a significantly better agreement with
results obtained by purely numerical methods. It can be clearly seen
that for certain classes of diagrams the asymptotic form of the vacuum
polarization function is not sufficient and power suppressed terms
play an important role and have to be included in the analysis.

\section*{Acknowledgements}
\label{sec:acknowledgements}
We like to thank K.G.~Chetyrkin and J.H.~K\"uhn for initiating the
project, fruitful discussions, and reading of the manuscript.  In
addition, we thank K.G.~Chetyrkin for valuable advice on computing the
vacuum polarization function in high-energy limit.  The work of
P.B. was supported by RFBR grant 11-02-01196. P.M. has been supported
in part by DFG Sonderforschungsbereich Transregio 9,
Computergest\"utzte Theoretische Teilchenphysik, and by the EU Network
{\sf LHCPHENOnet} PITN-GA-2010-264564.



%%%%%%%%%%%%%%%%%%%%%
%  Referenzen
%%%%%%%%%%%%%%%%%%%%%

\bibliography{biblio}{}
\bibliographystyle{elsarticle-num}
\end{document}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 
