\NeedsTeXFormat{LaTeX2e}
\documentclass[11pt,a4paper]{article}
\usepackage{jheppub}

\input{macros}
\allowdisplaybreaks
\newcommand{\ice}[1]{\relax}


\parindent 0pt
\sloppy
\parskip 6pt


\title{Top-Yukawa effects on the $\beta$-function of the strong coupling in the SM at four-loop level}
\author[a]{M. F. Zoller}
\affiliation[a]{Institut f\"ur Theoretische Teilchenphysik, Karlsruhe
  Institute of Technology (KIT), \mbox{D-76128 Karlsruhe, Germany}}
\emailAdd{max.zoller@kit.edu}

\abstract{We present analytical results for the QCD $\beta$-function extended to the gaugeless limit
of the unbroken phase of the Standard Model at four-loop level. Apart from the strong coupling itself
we include the top-Yukawa contribution and the Higgs self-coupling. We observe a non-naive
$\gaf$ contribution at order $\yt^4\gs^4$, a feature not encountered in lower loop orders.
}

\keywords{Renormalization Group, Standard Model, QCD}
\arxivnumber{}
\subheader{TTP15-031 }

\begin{document}
\maketitle

\section{Introduction}

An important feature of the perturbative treatment of any quantum field theory is the evolution of couplings, fields and masses with
the renormalization scale $\mu$, which is usually set to a characteristic energy scale of the physical process under consideration.
This evolution is described by the Renormalization Group (RG) functions, i.~e.~$\beta$-functions for the couplings and anomalous dimensions for fields and masses.

The $\beta$-function for any coupling $X$ is defined as
\be
\beta_{\sss{X}}(X,X_1,X_2,\ldots)=\mu^2\f{d X}{d \mu^2}    =\sum \limits_{n=1}^{\infty} \f{1}{(16\pi^2)^{n}}\,\beta_{\sss{X}}^{(n)}
\ee
and the anomalous dimension of a field $f$ as
\be
\gamma^{\sss{f}}_2=-\mu^2\f{d \text{ln} Z_{\sss{f}}^{-1}}{d \mu^2}    =
\sum \limits_{n=1}^{\infty} \f{1}{(16\pi^2)^{n}}\,\gamma_2^{{\sss{f}}\,(n)}
{},
\ee
where $Z_{\sss{f}}$ is the field strength renormalization constant\footnote{Note that we take the inverse of $Z_{\sss{f}}$ in the definition of
$\gamma^{\sss{f}}_2$. For an $n$-point vertex $V$ or a mass $m$ the anomalous dimension is defined as \mbox{$\gamma^{\sss{V}}_{n}=-\mu^2\f{d \text{ln} Z_{\sss{V}}}{d \mu^2}$} or
\mbox{$\gamma_{\sss{m}}=-\mu^2\f{d \text{ln} Z_{\sss{m}}}{d \mu^2}$}.} for the respective field. These functions are
power series in all couplings \mbox{$X,X_1,X_2,\ldots$} of the theory.
The $\beta$-functions and the anomalous dimensions of masses are independent of all gauge parameters $\xi$ whereas the anomalous dimensions of the fields are
not.
% 
% Note that the $\beta$-functions $\beta_{\sss{\gs}}$ and $\beta_{\sss{\yt}}$ are proportional to $\gs$ and $\yt$ respectively whereas
% $\beta_{\sss{\lambda}}$ has one part proportional to $\lambda$ and one part proportional to $\yt^4$ with no $\lambda$-dependence at all.

Recently the RG functions of the Standard Model (SM) were computed at three-loop accuracy. In the \msbar scheme $\beta$-functions do not depend on masses \cite{Collins:1974da}, hence
they can be computed in the unbroken phase of the SM.
For the gauge couplings $\gs,\gw$ and $\gb$ of the $\text{SU}_{\sss{C}}(3)$, $\text{SU}_{\sss{L}}(2)$ and
$\text{U}_{\sss{Y}}(1)$ subgroups of the SM the results were first published in \cite{PhysRevLett.108.151602,Mihaila:2012pz} and independently confirmed in \cite{Bednyakov:2012rb}.
For the top-Yukawa coupling $\yt$, which is the numerically largest Yukawa coupling by far, and the parameters of the Higgs potential 
$\lambda$ and $m^2$ the $\beta$-functions were first computed in the gaugeless limit, i.e. $\gw, \gb \to 0$, along with the anomalous dimensions of the fields involved \cite{Chetyrkin:2012rz}. 
Later $\beta_{\lambda}$ and $\beta_{m^2}$ were extended to the full SM \cite{Chetyrkin:2013wya}, confirmed by \cite{Bednyakov:2013eba,Bednyakov:2013cpa}, as well as
$\beta_{\yt}$ \cite{Bednyakov:2012en}, where the $\beta$-functions for the smaller Yukawa couplings were also added. 
The one- and two-loop $\beta$-functions for the gauge couplings
\cite{PhysRevLett.30.1343,PhysRevLett.30.1346,Jones1974531,Tarasov:1976ef,PhysRevLett.33.244,
Egorian:1978zx,PhysRevD.25.581,Fischler:1981is,Jack1985472,
Machacek198383,2loopbetayukawa,Ford:1992pn}, Yukawa couplings
\cite{Fischler1982385,Jack1985472,Machacek1984221,2loopbetayukawa} and Higgs potential parameters
\cite{Machacek198570,Jack1985472,2loopbetayukawa,Ford:1992pn}
have been known for a long time as well as partial three-loop results
\cite{Curtright:1979mg,Jones:1980fx,Tarasov:1980au,Tarasov1980429,3loopbetaqcd,Steinhauser:1998cm,Pickering:2001aq}. 
At four-loop level only the QCD $\beta$-function, i.e. $\beta_{\gs}(\gs)$ or equivalently 
$\beta_{\als}(\als)=2 \f{\als}{\gs}\beta_{\gs}$ with $\als=\f{\gs^2}{4\pi}$ is known \cite{4loopbetaqcd,Czakon:2004bu}.

Especially the evolution of the quartic Higgs self-coupling has received a lot of interest because of its close connection to the question of vacuum
stability in the Standard Model. It has been shown that the stability
of the SM vacuum up to some large energy scale $\Lambda \sim M_{\sss{Planck}}$ is approximately
equivalent to the requirement that the running coupling
$\lambda(\mu)>0$ for $\mu \leq \Lambda$ \cite{Cabibbo:1979ay,Ford:1992mv,Altarelli1994141}. 
The function $\beta_{\lambda}$ describing this evolution depends on all SM couplings an especially the large couplings $\yt$ and $\gs$ have a strong
influence. As the evolution of all couplings is interdependent a precision calculation for the evolution of all - at least of the five largest ($\gs$, $\yt$, $\gw$, $\gb$ and $\lambda$) -
is well motivated. Many analyses of this question have been performed
\cite{Bezrukov:2009db,Holthausen:2011aa,EliasMiro:2011aa,Xing:2011aa,Bezrukov:2012sa,Degrassi:2012ry,Chetyrkin:2012rz,Zoller:2012cv,Masina:2012tz,Zoller:2014cka,Zoller:2014xoa,Zoller:2013mra,Buttazzo:2013uya,Bednyakov:2015sca}
during the last years. 

In this paper we extend the QCD $\beta$-function to the gaugeless limit of the SM, i.~e.~we include the dependence on the top-Yukawa coupling $\yt$ and
the quartic Higgs self-coupling $\lambda$. This can be seen as a first step to all three gauge coupling $\beta$-functions in the full SM. 
To start with the gaugeless limit seems reasonable, first because at the energy-scales of our experiments $\yt$ is the second largest coupling in the SM after $\gs$,
followed by $\gw$, $\gb$ and $\lambda$. In order to renormalize fermion loops with four scalar legs we should also add counterterm $\propto \Phi^4$ to the Lagrangian of our simplified model.
This is exactly a contribution to the renormalization of $\lambda$ which makes it natural to include $\lambda$ as well. 

Secondly, the gaugeless limit of the SM provides an excellent opportunity to study the proper treatment of $\gaf$, which is introduced in the Yukawa-part of the Lagrangian.
This matrix is not well-defined in $D=4-2\eps$ dimensions and hence constitutes a non-trivial challenge.

The paper is structured as follows: In the following section the technical details, especially the treatment of $\gaf$, as well as
the automation of the calculation are discussed. Then the results are given and the relevance of the four-loop terms numerically determined at
the scale of the top quark mass. 

{\bf Note:} During the finishing process of this paper a similar calculation was published by another group  \cite{Bednyakov:2015ooa}.
Their calculation was not performed with massive tadpole integrals but rather with massless propagator-like integrals and in the Background field gauge. 
Both results achieved with different methods agree but for one term which is numerically small.

\section{Details of the calculation}

\subsection{Gaugeless limit of the SM}

The Lagrangian of the SM in the unbroken phase can be decomposed into
\be
\ssL_{\ssst{SM}}=\ssL_{\ssst{SU(3)$\times$SU(2)$\times$U(1)}}+\ssL_{\ssst{Yukawa}}+\ssL_{\sss{\Phi}}{}, \label{LSM}
\ee
where $\ssL_{\ssst{SU(3)$\times$SU(2)$\times$U(1)}}$ contains the kinetic terms of the fermions and gauge bosons,
their interactions and the necessary gauge fixing and ghost terms. The Yukawa part $\ssL_{\ssst{Yukawa}}$ describes
the coupling of the fermions to a scalar SU(2) doublet $\Phi=\vv{\Phi_1}{\Phi_2}$ which results in fermion masses and the coupling 
of fermions to the Higgs boson after Spontaneous Symmetry Breaking as well as the mixing of the quark generations. 
The scalar part $\ssL_{\sss{\Phi}}$
contains the kinetic term for the scalar field $\Phi$, its potential and its coupling to the electroweak gauge bosons
through the covariant derivative. 
In the gaugeless limit we neglect two smaller gauge couplings $\gw$ and $\gb$ (electroweak sector).
We also approximate the small Yukawa couplings, i.~e.~all but the top-Yukawa coupling $\yt$, by zero and arrive at a simplified model
which includes QCD and top-Yukawa effects as well as the scalar potential:
\be
\ssL=\ssL_{\sss{QCD}}+\ssL_{\sss{\yt}}+\ssL_{\sss{\Phi}} \label{Lgaugeless}
\ee
with
\bea
\ssL_{\sss{QCD}}&=&-\f{1}{4}G^a_{\mu \nu} G^{a\,\mu \nu}-\f{1}{2 (1-\xi)}\lb\p_\mu A^{a\,\mu}\rb^2 
+\p_\mu \bar{c}^a \p^{\mu}c^a+\gs f^{abc}\,\p_\mu \bar{c}^a A^{b\,\mu} c^c \nonumber \\
&+&\sum\limits_q\left\{\f{i}{2}\bar{q}\overleftrightarrow{\slashed{\p}}q+ \gs \bar{q}\slashed{A}^a T^a q\right\}{},
\label{LQCD} \\
\ssL_{\sss{\yt}}
&=&-\yt \left\{ \lb \bar{t} P_{\sss{R}} t\rb \Phi^{*}_2+\lb\bar{t} P_{\sss{L}} t\rb \Phi_2
-\lb\bar{b} P_{\sss{R}} t\rb \Phi^{*}_1-\lb \bar{t} P_{\sss{L}} b\rb \Phi_1
\right\} {},
\label{LYuk} \\
\ssL_{\sss{\Phi}}&=&\p_\mu \Phi^\dagger \p^\mu \Phi-m^2\Phi^\dagger\Phi-\lambda \lb\Phi^\dagger\Phi\rb^2{}. 
\label{LPhi} 
\eea  
Here $q$ runs over all quark flavours, the gluon field strength tensor is given by
\be
G^a_{\mu \nu}=\p_\mu A^a_\nu - \p_\nu A^a_\mu + \gs f^{abc}A^b_\mu A^c_\nu
\ee
and $f^{abc}$ are the structure constants of the colour gauge group with the generators $T^a$ which satisfy
\be \left[ T^a,T^b \right]=if^{abc}T^c.\ee

The Yukawa sector mixes left-handed (L) and right-handed (R) Weyl spinors which can be projected out from Dirac spinors
used in our Feynman rules by the application of the projectors
\be 
P_{\sss{L}}=\f{1}{2}\lb 1-\gamma_5\rb \qquad P_{\sss{R}}=\f{1}{2}\lb 1+\gamma_5\rb
{}.
\ee
The left- and right-handed parts of the quark fields and vertices participating in the Yukawa interaction are renormalized
differently.

The Lagrangian \eqref{Lgaugeless} is renormalized with the counterterms
\bea
\delta\!\ssL_{\sss{QCD}}&=&-\f{1}{4}\delta\! Z^{(2g)}_3 \lb \p_\mu A^a_\nu - \p_\nu A^a_\mu \rb^2
-\f{1}{2}\delta\! Z^{(3g)}_1 \gs f^{abc}\lb \p_\mu A^a_\nu - \p_\nu A^a_\mu \rb A^b_\mu A^c_\nu \nonumber\\
&-&\f{1}{4}\delta\! Z^{(4g)}_1 \gs^2 \lb f^{abc} A^b_\mu A^c_\nu \rb^2
+\delta\! Z^{(2c)}_3 \p_\mu \bar{c}^a \p^{\mu}c^a+\delta\! Z^{(ccg)}_1 \gs f^{abc}\,\p_\mu \bar{c}^a A^{b\,\mu} c^c \label{LQCDc}\\
&+&\sum\limits_q\left\{
\f{i}{2}\bar{q}\overleftrightarrow{\slashed{\p}}\left[\delta\! Z^{(2q)}_{2,L}P_{\sss{L}}+\delta\! Z^{(2q)}_{2,R}P_{\sss{R}}\right]q
+\gs \bar{q}\slashed{A}^a T^a\left[\delta\! Z^{(qqg)}_{1,L}P_{\sss{L}} +\delta\! Z^{(qqg)}_{1,R}P_{\sss{R}}\right]q
\right\}{}, \nonumber \\
\delta\!\ssL_{\sss{Yukawa}}&=&-\delta\! Z^{(tb\Phi)}_1 
\yt\left\{ \lb \bar{t} P_{\sss{R}} t\rb \Phi^{*}_2+\lb\bar{t} P_{\sss{L}} t\rb \Phi_2
-\lb\bar{b} P_{\sss{R}} t\rb \Phi^{*}_1-\lb \bar{t} P_{\sss{L}} b\rb \Phi_1
\right\} {},
\label{LYukc} \\
\delta\!\ssL_{\sss{\Phi}}&=&\delta\! Z_2^{(2\Phi)}\p_\mu \Phi^\dagger \p^\mu \Phi-m^2\, \delta\! Z_{\Phi^2}\Phi^\dagger\Phi
+\delta\! Z_1^{(4\Phi)}\lb\Phi^\dagger\Phi\rb^2{}.
\label{LPhic} 
\eea
All these renormalization constants were computed at three-loop level in the course of the calculations in \cite{Chetyrkin:2012rz}.
The simplest way to derive the renormalization constant for the strong gauge coupling $\gs$ is via
\be 
Z_{\gs}=\f{Z^{(ccg)}_{1}}{Z^{(2c)}_{3}\sqrt{Z^{(2g)}_{3}}} \label{Zgscomputation}
\ee
where we use the renormalization constants $Z=1+\delta Z$ in the \msbar-scheme. 
All divergent integrals are regularized in $D=4-2 \eps$ space time dimensions.

\subsection{Automation and calculation with massive tadpoles}

The calculation begins with the generation of all necessary 1PI Feynman diagrams with two external ghost or gluon legs
for $Z^{(2c)}_3$ or $Z^{(2g)}_3$ and with two external ghost and one external gluon leg for $Z^{(ccg)}_1$. This was done
with the program QGRAF \cite{QGRAF}. 

The C++ programs Q2E and EXP \cite{Seidensticker:1999bb,Harlander:1997zb} are then used to identify the topology of the diagram. Later we will Taylor expand in the external
momenta and use projectors on the integrals in order to make them scalar. For example the ghost-gluon vertex corrections are
proportional to the outgoing ghost momentum $q^\mu$, where $\mu$ is the Lorentz index of the gluon leg. Hence we expand
to first order in $q$, use the projector $\f{q^\mu}{q^2}$ on the integral and set $q \to 0$ after that.
This is allowed as \msbar{} renormalization constants do not depend on external momenta. After having set all external momenta to zero
we are left with tadpole integrals. 
The fermion traces, the expansion in the external momenta and the insertion of counterterms in one-loop, 
two-loop and three-loop diagrams was performed using FORM \cite{Vermaseren:2000nd,Tentyukov:2007mu}. 
The colour factors were computed with the \mbox{FORM} package COLOR \cite{COLOR}. 
The tadpole integrals up to three-loop order were computed with
the \mbox{FORM}-based package \mbox{MATAD}\cite{MATAD}.

At four-loop level there are two independent tadpole topologies, see Fig.~\ref{Topotad4l}.
\footnote{All Feynman diagrams in this paper have been drawn with the 
Latex package Axodraw \cite{Vermaseren:1994je}.}
\begin{figure}[!ht]
\begin{center}
  \begin{tabular}{ll}
    \begin{picture}(200,120) (0,0)
    \SetWidth{0.5}
    \SetColor{Black}
    \ArrowLine(0,0)(0,100) 
    \ArrowLine(0,100)(150,100)
    \ArrowLine(150,100)(150,0) 
    \ArrowLine(150,0)(0,0)
    \ArrowLine(50,0)(50,50)
    \ArrowLine(100,0)(100,50)
    \ArrowLine(50,50)(50,100)
    \ArrowLine(100,50)(100,100)
    \ArrowLine(50,50)(100,50)
    \Vertex(50,0){2}
    \Vertex(50,50){2}
    \Vertex(50,100){2}
    \Vertex(100,0){2}
    \Vertex(100,50){2}
    \Vertex(100,100){2}
    \Text(-5,50)[rc]{\Large{{$p_1$}}}  
    \Text(75,108)[cc]{\Large{{$p_2$}}}
    \Text(155,50)[lc]{\Large{{$p_3$}}}
    \Text(75,8)[cc]{\Large{{$p_4$}}}
    \Text(45,25)[rc]{\Large{{$p_5$}}}
    \Text(45,75)[rc]{\Large{{$p_6$}}}
    \Text(105,75)[lc]{\Large{{$p_7$}}}
    \Text(105,25)[lc]{\Large{{$p_8$}}}
    \Text(75,58)[cc]{\Large{{$p_9$}}}
  \end{picture} &
      \begin{picture}(200,120) (0,0)
    \SetWidth{0.5}
    \SetColor{Black}
    \ArrowLine(0,0)(0,100) 
    \ArrowLine(0,100)(150,100)
    \ArrowLine(150,100)(150,0) 
    \ArrowLine(150,0)(0,0)
    \ArrowLine(50,0)(50,50)
    \ArrowLine(100,0)(100,50)
    \Line(50,50)(74,74)
    \Line(100,50)(75,75)
    \ArrowLine(76,76)(100,100)
    \ArrowLine(75,75)(50,100)    
    \ArrowLine(50,50)(100,50)
    \Vertex(50,0){2}
    \Vertex(50,50){2}
    \Vertex(50,100){2}
    \Vertex(100,0){2}
    \Vertex(100,50){2}
    \Vertex(100,100){2}
    \Text(-5,50)[rc]{\Large{{$p_1$}}}  
    \Text(75,108)[cc]{\Large{{$p_2$}}}
    \Text(155,50)[lc]{\Large{{$p_3$}}}
    \Text(75,8)[cc]{\Large{{$p_4$}}}
    \Text(45,25)[rc]{\Large{{$p_5$}}}
    \Text(60,80)[rc]{\Large{{$p_6$}}}
    \Text(90,80)[lc]{\Large{{$p_7$}}}
    \Text(105,25)[lc]{\Large{{$p_8$}}}
    \Text(75,58)[cc]{\Large{{$p_{10}$}}}
  \end{picture} \\ planar topology: tad4lp & non-planar topology: tad4lnp
  \end{tabular} \end{center}
  \caption{Four-loop tadpole topologies: $p_1$, $p_2$, $p_3$, $p_4$ are independent loop momenta, the others are linear combinations
\mbox{$p_5=p_4-p_1$},
\mbox{$p_6=p_2-p_1$},
\mbox{$p_7=p_3-p_2$},
\mbox{$p_8=p_3-p_4$},
\mbox{$p_9=p_4-p_2$} and
\mbox{$p_{10}=p_4+p_2-p_1-p_3$}.}
  \label{Topotad4l}
\end{figure} 
All scalar products $p_i\cdotp_j$ ($i,j=1,\ldots,10$) can be written as linear combinations of the $p_i^2$
which can be expressed in terms of the scalar propagators $D_i=\f{1}{i}\f{1}{M^2-p_i^2}$ and the auxiliary
Mass $M^2$ (see below). Hence all four-loop integrals can be written in terms of functions
\be \text{TAD4l}(n_1,\ldots,n_{10}):=\int\!\mathrm{d}^D\!p_1\int\!\mathrm{d}^D\!p_2\int\!\mathrm{d}^D\!p_3\int\!\mathrm{d}^D\!p_4
\prod\limits_{i=1}^{10} D_i^{n_i}{}. \label{unredInt}\ee

The integrals \eqref{unredInt} can be reduced to Master Intgrals (MI) using 
FIRE \cite{Smirnov:2008iw}. For the huge number of integrals in such a calculation the C++ version of
FIRE 5 \cite{Smirnov:2014hma} is necessary. All MI needed for this computation can be found in \cite{Czakon:2004bu}.
The program FIESTA 3 \cite{Smirnov:2013eza} was used to numerically cross check these MI and some unreduced
integrals as a check for our setup.

In order to compute the divergent part of the needed self-energies and vertex corrections we use
the same method as in our previous calculations \cite{Chetyrkin:2012rz,Chetyrkin:2013wya}.
This method was suggested in \cite{Misiak:1994zw} and further developed in \cite{beta_den_comp}. 
A step-by-step explanation of this method can be found in \cite{Zoller:2014xoa}.
An auxiliary mass parameter $M^2$ is introduced in every propagator denominator.
A naive Taylor expansion in the external momenta is performed before applying the projector to scalar integrals. 
After that all external
momenta are set to zero which leaves us with scalar tadpole integrals. Subdivergences $\propto M^2$ are canceled by counterterms
\be \begin{split}
\f{M^2}{2}\delta\!Z_{\sss{M^2}}^{(2g)}\,A_\mu^a A^{a\,\mu}\;
% \f{M^2}{2}\delta\!Z_{\sss{M^2}}^{(2W)}\,W_\mu^a W^{a\,\mu},\; 
% \f{M^2}{2}\delta\!Z_{\sss{M^2}}^{(2B)}\,B_\mu B^{\mu}\;
\text{and} \;
\f{M^2}{2}\delta\!Z_{\sss{M^2}}^{(2\Phi)}\, \Phi^\dagger\Phi{}.
\end{split} \ee
which are computed order by order in perturbation theory and inserted in lower loop diagrams.
Note that this method only yields the divergent part of the integrals correctly, not the finite part.

\subsection{Treatment of $\gaf$}
The most important issue of this calculation is the proper treatment of $\gaf$ in dimensional regularization. In $D=4$ dimensions it can be defined as 
\be \gamma_5=i\gamma^0\gamma^1\gamma^2\gamma^3=\f{i}{4!}\eps_{\mu\nu\rho\sigma} \gamma^\mu\gamma^\nu\gamma^\rho\gamma^\sigma \text{ with }
 \eps_{0123}=1=-\eps^{0123} \label{gamma5} 
{}.
\ee
In most diagrams a naive treatment of $\gaf$ is sufficient, i.~e.~we use $\left\{\gaf,\gamma^\mu\right\}=0$ and $\gaf^2=1$,
valid in $D=4$ dimensions, until only
one or no $\gaf$ matrix remain on each fermion line, then discard diagrams with at least one $\gaf$.
This is valid for fermion lines with less than four Lorentz indices and momenta flowing into the fermion line.
\begin{figure}
 \begin{tabular}{cc}
 \scalebox{0.8}{
  \begin{picture}(230,45) (0,0) 
\ArrowArc(100,0)(50,180,360)
\ArrowArc(100,0)(50,0,180)
\Gluon(50,0)(0,0){6}{3}
\Gluon(150,0)(200,0){6}{3}
\Gluon(65,-35)(45,-75){6}{3}
\Gluon(135,-35)(165,-75){6}{3}
\DashArrowLine(45,75)(65,35){4}
\DashArrowLine(165,75)(135,35){4}
\SetColor{Red}
\Vertex(65,35){4}
\LongArrow(35,75)(55,35)
\LongArrow(175,75)(145,35)
 \Text(75,15)[cc]{\Large{\Red{$\gamma_5$}}}
\Text(0,-16)[cc]{{\Red{$\mu_1$}}}
\Text(200,-16)[cc]{{\Red{$\mu_2$}}}
\Text(45,-91)[cc]{{\Red{$\mu_3$}}}
\Text(165,-91)[cc]{{\Red{$\mu_4$}}}
\Text(25,85)[cc]{{\Red{$k_2$}}}
\Text(185,85)[cc]{{\Red{$k_1$}}}
\end{picture}}
& \scalebox{0.8}{
\begin{picture}(230,45) (0,0) 
 \ArrowLine(0,0)(200,0)
 \ArrowLine(180,0)(230,0)
 \ArrowLine(0,0)(50,0)
 \Gluon(50,0)(50,40){4}{4}
 \Gluon(80,0)(80,40){4}{4}
 \Gluon(120,0)(120,40){4}{4}
 \Gluon(180,0)(180,40){4}{4}
 \DashArrowLine(150,0)(150,40){4}
 \COval(115,50)(20,75)(0){Black}{Blue}
\Text(100,15)[cc]{{\Black{$\cdots$}}}
\SetColor{Red}
\Vertex(150,0){4}
\LongArrowArc(180,10)(30,220,320)
 \Text(180,-35)[cc]{\Large{\Red{$\gamma_5$}}}
\end{picture}}\vspace{20mm}\\
(a) & (b)
\end{tabular}
\caption{$\gaf$ on internal (a) and external (b) fermion lines}
\label{gammafiveintext}
\end{figure}
Fig.~\ref{gammafiveintext} shows the schematic cases of $\gaf$ appearing on internal and external fermion lines.
We start with the case of internal lines (see Fig.~\ref{gammafiveintext} (a)). In fact, for the calculation presented in this paper no external fermion lines appear.

As we set all momenta external to the whole Feynman diagram to zero for the computation of the UV divergent part
of the diagram external momenta to a fermion line ($k_1$, $k_2$,$\ldots$) are loop momenta from other loops.
Taking the trace over the closed fermion loop in $D=4$ dimensions yields a result with terms    
proportional to $\eps_{\mu_1\mu_2\mu_3\mu_4}$ and $\eps_{\mu_1 \mu_2 \alpha \beta}\, k_1^\alpha\, k_2^\beta$ and so on.
In order for the $\eps$-tensors not to vanish at least 4 free Lorentz structures are needed.
Else the diagram can be set to zero before we even employ dimensional regularization.

If we have only one internal fermion line with one $\gaf$ on it and the final result is known to be scalar (not pseudoscalar), as are the counterterms we want to compute here,
we can discard these terms as well. The only possibility for a non-naive contribution to the final result can appear in the case
of two (or more) fermion lines. Here the two $\eps$-tensors can be contracted and expressed in terms of the metric tensor
\be \eps^{\mu_1\mu_2\mu_3\mu_4} \eps_{\nu_1\nu_2\nu_3\nu_4}
%=-g^{[\mu_1}_{\,\,\,\,\nu_1} g^{\mu_2}_{\,\,\,\,\nu_2} g^{\mu_3}_{\,\,\,\,\nu_3} g^{\mu_4]}_{\,\,\,\,\nu_4} 
=-\sum\limits_{\pi}\,\text{sgn}(\pi) g^{\mu_{\pi(1)}}_{\,\,\,\,\nu_1} g^{\mu_{\pi(2)}}_{\,\,\,\,\nu_2} g^{\mu_{\pi(3)}}_{\,\,\,\,\nu_3} 
g^{\mu_{\pi(4)}}_{\,\,\,\,\nu_4}
{}, \label{epscontraction} \ee
where the sum is taken over all permutations $\pi$ of (1,2,3,4) and \be \text{sgn}(\pi)=\begin{cases}
                                                                                      +1 & \text{ for $\pi$ even} \\ -1 & \text{ for $\pi$ odd}
                                                                                     \end{cases}
{}. \ee
The lhs of \eqref{epscontraction} is composed of intrinsically four-dimensional objects whereas the rhs can be used in $D=4-2\eps$ dimensions,
introducing an uncertainty of $\mathcal{O}(\eps)$. However, if the integrals appearing in the calculation of the Feynman diagram in question have only
$\f{1}{\eps}$ poles the divergent part, which we are interested in here, is unaffected.

For completeness we want to make a short remark about external fermion lines, such as the one shown in Fig.~\ref{gammafiveintext} (b), as well.
Here we can anticommute the $\gaf$ to the end of the fermion line and hence outside of all loops. But if we use a projector on the external fermion line
in order to make the integral scalar and this involves taking a trace over the fermion line we have to treat it the same way as the
internal ones. In the case of the three-loop $\beta$-function for the Yukawa couplings a non-naive $\gaf$ effect from the contraction
of the $\eps$-tensors from an internal and an external fermion line was observed \cite{Chetyrkin:2012rz}.

In the the calculations needed for the renormalization constants in \eqref{Zgscomputation} only one type of diagram features two fermion lines with four external
Lorentz indices or loop-momenta to them, namely in the gluon propagator, when each external leg is attached to a different fermion loop and the two fermion loops are
connected by a gluon and two $\Phi$-lines. A planar example is shown in Fig.~\ref{diasggG5ep}

\begin{figure}[!ht]
  \begin{center}
    \begin{picture}(200,110) (0,-50)
    \SetWidth{0.5}
    \SetColor{Black}
    \Gluon(0,0)(30,0){4}{3}
    \ArrowLine(30,0)(70,40)
    \ArrowLine(70,-40)(30,0)
    \ArrowLine(70,40)(70,0)
    \ArrowLine(70,0)(70,-40)
        \Gluon(70,40)(120,40){4}{5}
        \DashArrowLine(70,0)(120,0){6}
        \DashArrowLine(120,-40)(70,-40){6}
    \ArrowLine(120,40)(160,0)
    \ArrowLine(160,0)(120,-40)
    \ArrowLine(120,0)(120,40)
    \ArrowLine(120,-40)(120,0)
    \Gluon(160,0)(190,0){4}{3}
\SetColor{Red}
        \Text(0,-10)[lt]{\LARGE{\Red{$\mu_1$}}}    
        \Text(190,-10)[rt]{\LARGE{\Red{$\mu_2$}}}
        \Vertex(30,0){3}
        \Vertex(160,0){3}
                \Text(70,45)[cb]{\LARGE{\Red{$\nu_1$}}}    
                \Text(120,45)[cb]{\LARGE{\Red{$\nu_2$}}}
                \Vertex(70,40){3}
                \Vertex(120,40){3}
\LongArrowArc(95,20)(12,210,160)
\Text(95,20)[cc]{\Large{\Red{$l_1$}}}
\LongArrowArc(95,-20)(12,210,160)
\Text(95,-20)[cc]{\Large{\Red{$l_2$}}}
\SetColor{Blue}
        \Vertex(70,0){3}
        \Vertex(120,0){3}
        \Vertex(70,-40){3}
        \Vertex(120,-40){3}
        \Text(65,0)[rc]{\Large{\Blue{$P_{\sss{L}}$}}}
        \Text(65,-45)[rc]{\Large{\Blue{$P_{\sss{R}}$}}}
        \Text(125,0)[lc]{\Large{\Blue{$P_{\sss{R}}$}}}
        \Text(125,-45)[lc]{\Large{\Blue{$P_{\sss{L}}$}}}
\SetColor{Black}
%     \Text(0,50)[lb]{\Large{\Black{(a)}}}
  \end{picture}
\end{center}
\caption{Diagram giving a non-naive $\gamma_5$-contribution to the gluon self-energy: each fermion line has two indices $\mu_i$ and $\nu_i$ (i=1,2) and two momenta
$l_1, l_2$, which can support a $\eps^{\mu_i \nu_i \alpha \beta}l_\alpha l_\beta$ term. The left- and right-handed projectors $P_{\sss{L,R}}$ introduce $\gaf$ into the diagram.} 
\label{diasggG5ep}
\end{figure}
It was checked explicitly that only $\f{1}{\eps}$ poles appear in the integrals needed for these diagrams. In fact, as an additional precaution we checked that
at $\mathcal{O}(\eps)$ completely antisymmetric and completely symmetric structures composed of the metric and the eight indices appearing in the $\eps$ tensors do not give contributions to
the divergent part. This was implemented as
\bea 
\eps^{\mu_1\mu_2\mu_3\mu_4} \eps_{\nu_1\nu_2\nu_3\nu_4} = 
&-&\sum\limits_{\pi}\,\text{sgn}(\pi) g^{\mu_{\pi(1)}}_{\,\,\,\,\nu_1} g^{\mu_{\pi(2)}}_{\,\,\,\,\nu_2} g^{\mu_{\pi(3)}}_{\,\,\,\,\nu_3} 
g^{\mu_{\pi(4)}}_{\,\,\,\,\nu_4}
\lb 1+\eps \cdot \text{C}^{\text{as}}\rb \nonumber\\
&+& \eps \cdot \text{C}^{\text{s}}\sum\limits_{\pi}\, g^{\mu_{\pi(1)}}_{\,\,\,\,\nu_1} g^{\mu_{\pi(2)}}_{\,\,\,\,\nu_2} g^{\mu_{\pi(3)}}_{\,\,\,\,\nu_3} g^{\mu_{\pi(4)}}_{\,\,\,\,\nu_4}{},
\eea
where the labels $\text{C}^{\text{as,s}}$ parametrize the uncertainty introduced through \eqref{epscontraction} being applied in $D=4-2\eps$. As they drop out in the divergent term of our final result
we are satisfied that $\gaf$ can be treated in this way. An alternative that was also checked is to use \eqref{epscontraction}
and an $\tilde{\eps}$ instead of $\eps$ in the fermion traces and check that $\f{1}{\eps^i}$ with $i>1$ drop out in the sum of all diagrams.
After that we set $\tilde{\eps}=\eps$.

\section{Results \label{res:beta}}
In this section we give the results for the four-loop $\beta$-function of the strong coupling $\gs$ in the gaugeless limit of the SM.
% All results of this work can found at\\
% \texttt{\bf http://www-ttp.particle.uni-karlsruhe.de/Progdata/ttp15/ttp15-.../}.
For a gerneric SU($\dF$) gauge group the colour factors are expressed
through the quadratic Casimir operators $\cf$ and $\ca$ of the 
fundamental and the adjoint representation of the corresponding Lie algebra.
The dimension of the fundamental representation is called $\dF$. The adjoint representation has dimension $\Ng$ and the trace
$\tr$ defined by \mbox{$\tr \delta^{ab}=\textbf{Tr}\lb T^a T^b\rb$}  
with the group generators $T^a$ of the fundamental representation. In addition we need a few higher order invariants
constructed from the symmetric tensors
\bea
d_{\scriptscriptstyle{F}}^{abcd} &=& \f{1}{6} \text{Tr} 
\lb T^a T^b T^c T^d + T^a T^b T^d T^c + T^a T^c T^b T^d \right. \nonumber \\
&+& \left. T^a T^c T^d T^b + T^a T^d T^b T^c + T^a T^d T^c T^b \rb{}.
\eea
from the generators of the fundamental representation and analogously $d_{\scriptscriptstyle{A}}^{abcd}$ constructed from the generators of the adjoint representation.
The combinations needed and their SU($\dF$) values are
\bea
\dFFfNg &=& \frac{\dF^4-6\dF^2+18}{96\dF^2},\;\;\;\;
\dFAfNg = \frac{\dF(\dF^2+6)}{48},\\
\dAAfNg &=& \frac{\dF^2(\dF^2+36)}{24}.\nonumber
\eea
Furthermore for SU($\dF$) we have
\be \tr = \f{1}{2},\;\;\;\; \cf = \f{\dF^2-1}{2\dF},\;\;\;\; \ca = \dF,\;\;\;\; \Ng = \dF^2-1{}.\ee
The number of active fermion flavours is denoted by $\Nf$ (=6 in the SM).
\be
\begin{split}
\f{\beta_{\sss{\gs}}^{(4)}}{\gs}= &
         \gs^8   \left(
           \f{40}{9} \dAAfNg
          - \f{150653}{972} \ca^4
          - \f{256}{9} \Nf \dFAfNg
          - 23 \Nf \tr \cf^3     \right. \\ & \left.
          + \f{2102}{27} \Nf \ca \tr \cf^2
          - \f{7073}{486} \Nf \ca^2 \tr \cf
          + \f{39143}{162} \Nf \ca^3 \tr
          + \f{352}{9} \Nf^2 \dFFfNg     \right. \\ & \left.
          - \f{676}{27} \Nf^2 \tr^2 \cf^2
          - \f{8576}{243} \Nf^2 \ca \tr^2 \cf
          - \f{3965}{81} \Nf^2 \ca^2 \tr^2
          - \f{616}{243} \Nf^3 \tr^3 \cf     \right. \\ & \left.
          - \f{212}{243} \Nf^3 \ca \tr^3
          - \f{352}{3} \zeta_{3} \dAAfNg
          + \f{22}{9} \zeta_{3} \ca^4
          + \f{832}{3} \zeta_{3} \Nf \dFAfNg     \right. \\ & \left.
          - \f{176}{9} \zeta_{3} \Nf \ca \tr \cf^2
          + \f{328}{9} \zeta_{3} \Nf \ca^2 \tr \cf
          - \f{68}{3} \zeta_{3} \Nf \ca^3 \tr
          - \f{256}{3} \zeta_{3} \Nf^2 \dFFfNg     \right. \\ & \left.
          + \f{352}{9} \zeta_{3} \Nf^2 \tr^2 \cf^2
          - \f{224}{9} \zeta_{3} \Nf^2 \ca \tr^2 \cf
          - \f{112}{9} \zeta_{3} \Nf^2 \ca^2 \tr^2
          \right) \\ &
           + \gs^6 \yt^2   \left(
          - 3 \tr \cf^2
          - \f{523}{18} \ca \tr \cf
          - \f{985}{9} \ca^2 \tr
          + \f{322}{9} \Nf \tr^2 \cf     \right. \\ & \left.
          + \f{218}{9} \Nf \ca \tr^2
          + 72 \zeta_{3} \tr \cf^2
          + 36 \zeta_{3} \ca \tr \cf
          \right) \\ &
                 + \gs^4 \yt^4   \left(
          - 3 \tr \cf
          + \f{41}{2} \tr \cf \dR
          + 36 \ca \tr
          + 25 \ca \tr \dR \right. \\ & \left.
          - 24 \zeta_{3} \tr \cf \dR    
          +{\bl \tr^2\lb
           \f{100}{3} 
          - 16 \zeta_{3}  \rb}
          \right) \\ &
                + \gs^2 \yt^6   \left(
          - \f{21}{4} \tr
          - 29 \tr \dR
          - \f{3}{2} \tr \dR^2
          - 6 \zeta_{3} \tr
          \right)  \\ &
       - 30 \gs^2 \yt^4 \lambda \tr  
        + 36 \tr \gs^2 \yt^2 \lambda^2 .
\end{split}
\label{4lbetags}
\ee
This is in agreement with \cite{Bednyakov:2015ooa} except for the term $\propto\gs^4 \yt^4 \tr^2$ (blue) which is the only one affected
by non-naive $\gaf$ contributions as explained above. The naive and non-naive (i.~e.~stemming from the contraction of two $\eps$-tensors)
contributions are
\be 
\gs^4 \yt^4  \tr^2 \lb \f{100}{3} 
          - 16 \zeta_{3} \rb
=  \gs^4 \yt^4  \tr^2  \lb     
          \underbrace{ 24 }_{\text{ (naive)}}
          + \underbrace{ \f{28}{3} }_{\text{ (non-naive)}}
          - \underbrace{ 48 \zeta_{3}}_{\text{ (naive)}}
          + \underbrace{ 32 \zeta_{3}}_{\text{ (non-naive)}} \rb. \ee
The lower loop results are
\bea
\f{\beta_{\sss{\gs}}^{(3)}}{\gs} &=&
           \gs^6   \left(
          - \f{2857}{108} \ca^3
          - \Nf \tr \cf^2
          + \f{205}{18} \Nf \ca \tr \cf  \right. \nonumber \\ & &\left.
          + \f{1415}{54} \Nf \ca^2 \tr
          - \f{22}{9} \Nf^2 \tr^2 \cf
          - \f{79}{27} \Nf^2 \ca \tr^2
          \right)\\ & &
           -  \gs^4 \yt^2   \left(
           3 \tr \cf
          + 12 \ca \tr
          \right) 
                +  \gs^2 \yt^4   \left(
          + \f{9}{2} \tr
          + \f{7}{2} \tr \dR
          \right), \nonumber \\
\f{\beta_{\sss{\gs}}^{(2)}}{\gs} &=&  
        \gs^4   \left(
          - \f{17}{3} \ca^2
          + 2 \Nf \tr \cf
          + \f{10}{3} \Nf \ca \tr
          \right)       -2 \gs^2 \yt^2, \\
\f{\beta_{\sss{\gs}}^{(1)}}{\gs} &=& 
        \gs^2   \left(
          - \f{11}{6} \ca
          + \f{2}{3} \Nf \tr
          \right).
\eea
in agreement with \cite{Chetyrkin:2012rz}.
The pure QCD part of \eqref{4lbetags} agrees with \cite{4loopbetaqcd,Czakon:2004bu}.

For convenience we also give the $\beta$-function for $\als$. We absorb the loop factor $\f{1}{16\pi^2}$ into
\be \as=\f{\gs^2}{(4\pi)^2}=\f{\als}{4\pi},\quad \at=\f{\yt^2}{(4\pi)^2},\quad \allam=\f{\lambda}{(4\pi)^2}{} \ee and define
\be \beta_{\als}(\as,\at,\allam)=\sum \limits_{n=1}^{\infty} \beta_{\sss{\als}}^{(n)}(\as,\at,\allam).\ee
We find
\be
\begin{split}
\f{\beta_{\sss{\als}}^{(4)}}{\als}=& 
        \as^4   \left(
           \f{80}{9} \dAAfNg
          - \f{150653}{486} \ca^4
          - \f{512}{9} \Nf \dFAfNg
          - 46 \Nf \tr \cf^3      \right. \\ & \left.
          + \f{4204}{27} \Nf \ca \tr \cf^2 
          - \f{7073}{243} \Nf \ca^2 \tr \cf
          + \f{39143}{81} \Nf \ca^3 \tr  \right. \\ & \left.
          + \f{704}{9} \Nf^2 \dFFfNg     
          - \f{1352}{27} \Nf^2 \tr^2 \cf^2
          - \f{17152}{243} \Nf^2 \ca \tr^2 \cf \right. \\ & \left.
          - \f{7930}{81} \Nf^2 \ca^2 \tr^2
          - \f{1232}{243} \Nf^3 \tr^3 \cf     
          - \f{424}{243} \Nf^3 \ca \tr^3
          - \f{704}{3} \zeta_{3} \dAAfNg \right. \\ & \left.
          + \f{44}{9} \zeta_{3} \ca^4 
          + \f{1664}{3} \zeta_{3} \Nf \dFAfNg     
          - \f{352}{9} \zeta_{3} \Nf \ca \tr \cf^2\cf\right. \\ & \left.
          + \f{656}{9} \zeta_{3} \Nf \ca^2 \tr 
          - \f{136}{3} \zeta_{3} \Nf \ca^3 \tr 
          - \f{512}{3} \zeta_{3} \Nf^2 \dFFfNg   \cf\right. \\ & \left.  
          + \f{704}{9} \zeta_{3} \Nf^2 \tr^2 \cf^2
          - \f{448}{9} \zeta_{3} \Nf^2 \ca \tr^2 \cf
          - \f{224}{9} \zeta_{3} \Nf^2 \ca^2 \tr^2
          \right) \\ &
       + \at \as^3   \left(
          - 6 \tr \cf^2
          - \f{523}{9} \ca \tr \cf
          - \f{1970}{9} \ca^2 \tr
          + \f{644}{9} \Nf \tr^2 \cf \right. \\ & \left.
          + \f{436}{9} \Nf \ca \tr^2
          + 144 \zeta_{3} \tr \cf^2
          + 72 \zeta_{3} \ca \tr \cf
          \right) \\ &
       + \at^2 \as^2   \left(
          - 6 \tr \cf
          + 41 \tr \cf \dR
          + 72 \ca \tr
          + 50 \ca \tr \dR  \right. \\ & \left.
          - 48 \zeta_{3} \tr \cf \dR
          + {\bl \tr^2 \lb  \f{200}{3} 
          - 32 \zeta_{3}  \rb}
          \right) \\ &
       + \at^3 \as   \left(
          - \f{21}{2} \tr
          - 58 \tr \dR
          - 3 \tr \dR^2
          - 12 \zeta_{3} \tr
          \right)\\ &
                 + \at^2 \as \allam   \left(
          - 60 \tr
          \right)        + \at \as \allam^2   \left(
          + 72 \tr
          \right),
\label{4lbetaas}
\end{split}
\ee
where
\be \at^2 \as^2 \tr^2 \lb + \f{200}{3} \tr^2 - 32 \zeta_{3} \tr^2 \rb
= \at^2 \as^2 \tr^2   \lb       \underbrace{ 48 }_{\text{(naive)}}
          + \underbrace{ \f{56}{3}  }_{     \text{(non-naive)}}
          - \underbrace{  96 \zeta_{3} }_{  \text{(naive)}}
          + \underbrace{  64 \zeta_{3}  }_{ \text{(non-naive)}} \rb. \ee
and
\bea
\f{\beta_{\sss{\als}}^{(3)}}{\als} &=& 
       + \as^3   \left(
          - \f{2857}{54} \ca^3
          - 2 \Nf \tr \cf^2
          + \f{205}{9} \Nf \ca \tr \cf \right. \nonumber \\ & &\left.
          + \f{1415}{27} \Nf \ca^2 \tr
          - \f{44}{9} \Nf^2 \tr^2 \cf
          - \f{158}{27} \Nf^2 \ca \tr^2
          \right)\\ & &
                 + \at \as^2   \left(
          - 6 \tr \cf
          - 24 \ca \tr
          \right)
                 + \at^2 \as   \left(
          + 9 \tr
          + 7 \tr \dR
          \right), \nonumber \\
\f{\beta_{\sss{\als}}^{(2)}}{\als} &=& 
       + \as^2   \left(
          - \f{34}{3} \ca^2
          + 4 \Nf \tr \cf
          + \f{20}{3} \Nf \ca \tr
          \right) 
                 + \at \as   \left(
          - 4 \tr
          \right), \\
\f{\beta_{\sss{\als}}^{(1)}}{\als} &=& 
       + \as   \left(
          - \f{11}{3} \ca
          + \f{4}{3} \Nf \tr
          \right) {}.
\eea

Now we want to give a numerical evaluation of the $\beta$-functions at the scale of the top mass in order to get an idea
of the size of the new terms.
For \mbox{$M_t \approx 173.34\pm 0.76$ GeV} \cite{ATLAS:2014wva}, \mbox{$M_H \approx 125.09\pm 0.24$ GeV}\cite{Aad:2015zhl} 
and \mbox{$\als(M_Z)=0.1184\pm 0.0007$} \cite{Bethke:2012jm}
we get the couplings in the $\overline{\text{MS}}$-scheme at this scale using two-loop matching relations \cite{Buttazzo:2013uya}
\bea 
\gs(M_t)&=&1.1666 \pm 0.0035 \text{(exp)}, \nonumber \\
\yt(M_t)&=&0.9369 \pm 0.0046 \text{(exp)} \pm 0.0005 \text{(theo)}, \\
\lambda(M_t)&=&0.1259 \pm 0.0005 \text{(exp)} \pm 0.0003 \text{(theo)} \nonumber
\eea
where the experimental uncertainty (exp) stems from $M_t, M_H$ and $\als(M_Z)$ and the theoretical one (theo) from the matching of on-shell to \msbar{} parameters
(these are taken from \cite{Buttazzo:2013uya}). We find\footnote{The labels under the braces indicate from which part
of the $\beta$-function the contributions come.}
\bea
\f{\beta_{\sss{\gs}}^{(2)}}{\beta_{\sss{\gs}}^{(1)}}
&=& \underbrace{3.20\times 10^{-2}}_{\gs^4} 
 \underbrace{+1.59\times 10^{-3}}_{\gs^2\yt^2}, \\
\f{\beta_{\sss{\gs}}^{(3)}}{\beta_{\sss{\gs}}^{(1)}}
&=& \underbrace{ -3.45\times 10^{-4} }_{\gs^6} 
\underbrace{ +2.74\times 10^{-4} }_{\gs^4\yt^2} 
\underbrace{ -6.62\times 10^{-5} }_{\gs^2\yt^4} , \\
\f{\beta_{\sss{\gs}}^{(4)}}{\beta_{\sss{\gs}}^{(1)}}
&=&
\underbrace{2.26  \times 10^{-4}}_{\gs^8}
\underbrace{+2.47  \times 10^{-5}}_{\gs^6\yt^2}
\underbrace{-1.06  \times 10^{-5}}_{\gs^4 \yt^4 (\text{naive})}
\underbrace{-9.09 \times 10^{-7}}_{\gs^4\yt^4 (\text{non-naive})}\\ & &
\underbrace{ +2.77 \times 10^{-6}}_{\gs^2\yt^6} 
\underbrace{+ 1.06 \times 10^{-7}}_{\gs^2\yt^4\lambda}
\underbrace{-1.82  \times 10^{-8}}_{\gs^2\yt^2\lambda^2}
\eea
We see that the top-Yukawa contributions have a sizable impact on the four-loop $\beta$-function for the strong coupling.
The part $\propto\gs^6\yt^2$ increases it by $\sim 10 \%$ and the part $\propto\gs^4\yt^4$ decreases it by 
$\sim 5 \%$ at this scale compared to the pure QCD contribution $\propto\gs^8$. The non-naive term
gives a $\sim 0.4 \%$ contribution. Numerically the discrepancy with \cite{Bednyakov:2015ooa} is therefore 
negligible in any phenomenological analysis. At some point it should be resolved nevertheless.

\section{Conclusions \label{last}}

We have presented an analytical result for the four-loop $\beta$-function of the strong coupling $\gs$ in the gaugeless limit of the SM.
This constitutes an important extension of the well-known QCD result as top-Yukawa coupling is numerically the next important coupling after $\gs$, 
at least at the electroweak scale. Furthermore, this is an important step towards a complete calculation of the four-loop $\beta$-functions of the
gauge couplings in the full SM.

An important feature of this result is the non-naive $\gaf$ contribution $\propto \gs^4 \yt^4$. 
In the pure gauge boson and fermion sector of the SM, given by $\ssL_{\ssst{SU(3)$\times$SU(2)$\times$U(1)}}$, all non-naive contributions cancel in the sum of all diagrams,
making this part of the SM anomaly free. This has been explicitly checked during the calculation of the three-loop $\beta$-functions for the gauge couplings in the SM \cite{PhysRevLett.108.151602,Mihaila:2012pz}.
Here we see that with the inclusion of a scalar field non-naive contributions may appear in higher orders and special care will have to be taken when a attempting a complete calculation 
of four-loop $\beta$-functions in the SM.

\section*{Acknowledgements}
I am grateful to K.~G.~Chetyrkin for many useful discussions during the course of this project and
helpful comments on this paper. 

I also thank  J.~H.~K\"uhn for his support. 

This work has been supported in part by the Deutsche Forschungsgemeinschaft in the
Sonderforschungsbereich/Transregio SFB/TR-9 ``Computational Particle
Physics''.


\bibliographystyle{JHEP}

\bibliography{LiteraturSM}



\end{document}








