format: text

text:
# CERN-TH/2002-379 DESY 02-223 hep-ph/0212294 December 2002 rhad: a program for the evaluation of the hadronic 𝑅-ratio in the perturbative regime of QCD

## 1 General structure of $`𝑹\mathbf{(}𝒔\mathbf{)}`$

We are considering the fully inclusive production of quark pairs in $`e^+e^{}`$ annihilation (for a review see Ref. ). The tree-level diagram for this process in shown in Fig. 1$`(a)`$. At leading order (LO) perturbative QCD (pQCD), the cross section as a function of the squared center-of-mass (c.m.s.) energy, $`s`$, has thresholds at the points $`s=4m_Q^2`$ with $`Q=d,u,s,c,b,t`$. However, close to threshold, fixed-order pQCD is no longer applicable. Below threshold, non-perturbative effects lead to the formation of bound states of the quark–anti-quark pair, which appear as more or less sharp peaks in the cross section.111Due to the large width no bound state is formed in the top quark system. However, a peak in the cross section remains around $`s=4m_t^2`$. In our perturbative framework of charm and bottom production, these narrow resonances cannot be described, and are therefore not included in rhad. In addition, we have to spare out the region between the physical threshold $`s_Q^{\mathrm{low}}`$ and the beginning of the more or less flat continuum region at $`s_Q^{\mathrm{thr}}`$, where the cross section exhibits rapid variations. In the case of charm production, for example, the limits would be $`\sqrt{s_c^{\mathrm{low}}}2m_D3.73`$ GeV and $`\sqrt{s_c^{\mathrm{thr}}}4.8`$ GeV, while for bottom production they are $`\sqrt{s_b^{\mathrm{low}}}2m_B10.52`$ GeV and $`\sqrt{s_b^{\mathrm{thr}}}11.2`$ GeV. Even though rhad provides a numerical value in these threshold regions, one should not try to give this value any physical significance; rhad will print a warning that reminds the user of this fact.

Furthermore, we have to exclude the low-energy region, below about $`\sqrt{s_{\mathrm{min}}}=1.8`$ GeV, where the validity of perturbation theory is doubtful. Thus, mass effects from $`u`$-, $`d`$\- and $`s`$-quarks are negligible, which is why we will consider these quarks as massless throughout the paper.

Instead of the cross section $`\sigma (e^+e^{}\mathrm{hadrons})`$, we will express the results in terms of the so-called hadronic $`R`$-ratio,

$$\begin{array}{c}\hfill R(s)=\frac{\sigma (e^+e^{}\mathrm{hadrons})}{\sigma (e^+e^{}\mu ^+\mu ^{})},\end{array}$$
(1)
where

$$\begin{array}{c}\hfill \sigma (e^+e^{}\mu ^+\mu ^{})=\frac{4\alpha ^2\pi }{3s}\end{array}$$
(2)
is the high-energy limit for the photon-mediated lowest order muon pair production. In this paper we concentrate on the contributions induced by photon exchange. $`Z`$-boson exchange is only relevant at high energies. There, however, apart from the QCD corrections222For a discussion of corrections to the axial-vector correlator see Ref. ., also the electro-weak corrections become important ; these will not be addressed in this paper. Nevertheless, rhad can be used to evaluate the vector contribution to the production of top quarks in the same way as for charm and bottom production.

We write $`R(s)`$ as a sum of contributions from individual quark flavors, plus a so-called singlet term:

$$\begin{array}{c}\hfill R(s)=\underset{Q=d,u,s,c,b,t}{}R_Q(s)+R_{\mathrm{sing}}(s).\end{array}$$
(3)
The “non-singlet” contributions $`R_Q(s)`$ are all proportional to the square of the respective quark charge. The (numerically small) singlet piece $`R_{\mathrm{sing}}`$ appears for the first time at order $`\alpha _s^3`$ and will be discussed in Sect. 6.

$`R_Q(s)`$ is defined to be zero below the continuum region of $`Q\overline{Q}`$ production ($`s_Q^{\mathrm{thr}}`$), i.e., we do not attempt to describe the resonance regime:

$$\begin{array}{c}\hfill R_Q(s<s_Q^{\mathrm{thr}})=0.\end{array}$$
(4)
Our main concern are the radiative corrections to $`R_Q(s)`$, predominantly arising from gluons exchanged between, or emitted from, the $`Q\overline{Q}`$ pair, but the lowest order QED effects, though very small, will also be included. Thus we write

$$\begin{array}{cc}\hfill R_Q(s)& =\theta (ss_Q^{\mathrm{thr}})\left[\underset{n0}{}\left(\frac{\alpha _s}{\pi }\right)^nR_Q^{(n)}(s)+\delta R_Q^{\mathrm{QED}}(s)\right].\hfill \end{array}$$
(5)
This formula generally describes the QCD/QED effects to $`R_Q(s)`$, if $`\delta R_Q^{\mathrm{QED}}(s)`$ contains all contributions that vanish when electro-magnetic corrections are ignored. Both $`\alpha _s`$ and $`R_Q^{(n)}(s)`$ explicitely depend on the renormalization scale $`\mu `$, while $`R_Q(s)`$ is invariant under $`\mu `$–variation, up to and including the order of perturbation theory that is being considered.

Since the QCD corrections only affect the quarks in Fig. 1$`(a)`$, the leptonic part of the diagrams can be factored out and one remains with corrections to the $`\gamma Q\overline{Q}`$ vertex, and real radiation of quarks and gluons. Beyond order $`\alpha _s`$, the most convenient way to evaluate the fully inclusive rate is to compute the photon self energy and relate it to the rate $`\gamma ^{}Q\overline{Q}+X`$ through the optical theorem, i.e., by taking the imaginary part in the time-like region of the external momentum. Explicitely,

$$\begin{array}{c}\hfill R(s)=12\pi \mathrm{Im}\mathrm{\Pi }(s+iϵ),\end{array}$$
(6)
where

$$\begin{array}{c}\hfill \mathrm{\Pi }(q^2)=\frac{1}{3q^2}\left(g_{\mu \nu }+\frac{q_\mu q_\nu }{q^2}\right)\mathrm{\Pi }_{\mu \nu }(q)\end{array}$$
(7)
is the transversal part of the photon polarization function, $`\mathrm{\Pi }_{\mu \nu }(q)`$.

The different contributions to $`R(s)`$, as we classify them in our paper, will be illustrated by sample diagrams contributing to $`\mathrm{\Pi }(q^2)`$ in what follows. The LO terms are thus represented by Fig. 1$`(b)`$. Since the imaginary part of these diagrams is obtained by the sum of all possible cuts, it includes real and virtual contributions simultaneously.

A final remark concerning the counting of loops is in order. We understand by the $`n`$-loop QCD contribution to $`R(s)`$ the sum of real and virtual corrections to order $`\alpha _s^n`$. Computing $`R(s)`$ via the optical theorem at $`n`$-loop level thus requires the evaluation of the imaginary part of the photon self-energy to $`n+1`$ loops.

## 2 Running and decoupling of $`\alpha _s`$ and the quark masses

As already mentioned in the previous section, the quantity $`R(s)`$ depends on the strong coupling constant $`\alpha _s`$, the heavy quark masses $`m_Q`$, and the renormalization scale $`\mu `$. All of them are specified in the subroutine parameters which is discussed in Appendix C. The dependence of $`R(s)`$ on $`\mu `$ is explicit in terms of logarithms, but also implicit as, for example, in the argument of $`\alpha _s`$. At higher orders in perturbation theory, one needs to specify the scheme in which $`\alpha _s`$ and the quark masses $`m_Q`$ are evaluated. For $`\alpha _s`$, we will always adopt the $`\overline{\text{MS}}`$ scheme which is commonly used in QCD calculations. Concerning the quark masses, it is generally more appropriate to use the pole mass in the energy region close to the production threshold of the quarks. On the other hand, the use of the running mass and setting $`\mu ^2=s`$ resums part of the potentially large logarithms $`\mathrm{ln}m_Q^2/s`$ in the high-energy region. For this reason, rhad provides the logical parameter lmsbar that allows the user to switch between the pole and the $`\overline{\text{MS}}`$ mass definition for the evaluation of $`R(s)`$. If not stated otherwise, we denote by $`m_Q`$ a generic quark mass in what follows. In those cases where it is necessary to distinguish between the pole and $`\overline{\text{MS}}`$ quark mass, we denote them by $`M_Q`$ and $`\overline{m}_Q\overline{m}_Q^{(n_f)}(\mu )`$, respectively. The conversion formula between $`M_Q`$ and $`\overline{m}_Q`$ is given in App. E.4, Eq. (67). The “scale-invariant” mass is defined recursively as $`\overline{m}_Q(\overline{m}_Q)\overline{m}_Q^{(n_f)}(\overline{m}_Q^{(n_f)})`$, where the number of active flavors is equal to $`n_f=4`$ for charm, $`n_f=5`$ for bottom, and $`n_f=6`$ for top quarks.

If one chooses to evaluate $`R(s)`$ in terms of $`\overline{\text{MS}}`$ masses, the input required by rhad is $`\alpha _s^{(5)}(M_Z)`$ and the scale invariant masses $`\overline{m}_Q(\overline{m}_Q)`$ ($`Q=c,b,t`$). In addition, one needs to specify the renormalization scale $`\mu `$ and the c.m.s. energy squared, $`s`$. From that, rhad will determine the number of active flavors $`n_f`$, as well as the parameters $`\alpha _s^{(n_f)}(\mu )`$ and $`\overline{m}_Q^{(n_f)}(\mu )`$, which will be used for the evaluation of $`R(s)`$. As already mentioned before, a natural choice (and the default value in rhad) for the renormalization scale is $`\mu ^2=s`$.

The number of active flavors is determined through the threshold variables $`s_c^{\mathrm{thr}}`$, $`s_b^{\mathrm{thr}}`$, and $`s_t^{\mathrm{thr}}`$. When $`s<s_c^{\mathrm{thr}}`$, it assumes the smallest possible value, $`n_f=3`$, and it increases by one every time one of the threshold variables is crossed, up to $`n_f=6`$ for $`s>s_t^{\mathrm{thr}}`$.

On the other hand, in the definition of $`\alpha _s^{(n_f)}(\mu )`$ one has the freedom to choose “matching scales” $`\mu _c`$, $`\mu _b`$, and $`\mu _t`$, at which the transition from $`n_f`$ to $`n_f\pm 1`$ is performed. For example, assume that $`n_f=4`$, and the renormalization scale is set to some specific value $`\mu =\overline{\mu }`$. The procedure to compute $`\alpha _s^{(4)}(\overline{\mu })`$ at $`n`$-loop order as implemented in rhad is as follows: In a first step the renormalization group equation (RGE) for the strong coupling (here and in what follows, we refer to App. E.4 for the definition of the coefficients $`\beta `$, $`\zeta _g`$, $`\gamma _m`$, and $`\zeta _m`$),

$`\mu ^2{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\mu ^2}}{\displaystyle \frac{\alpha _s^{(n_f)}(\mu )}{\pi }}`$ $`=`$ $`\beta ^{(n_f)}\left(\alpha _s^{(n_f)}(\mu )\right)={\displaystyle \underset{i=0}{\overset{n1}{}}}\beta _i^{(n_f)}\left({\displaystyle \frac{\alpha _s^{(n_f)}(\mu )}{\pi }}\right)^{i+2},`$ (8)
is solved numerically in order to obtain $`\alpha _s^{(5)}(\mu _b)`$ from the input value $`\alpha _s^{(5)}(M_Z)`$. The use of the decoupling relation

$`\alpha _s^{(n_f1)}(\mu _b)`$ $`=`$ $`\left(\zeta _g\right)^2\alpha _s^{(n_f)}(\mu _b)`$ (9)
to $`(n1)`$-loop order, leads to $`\alpha _s^{(4)}(\mu _b)`$. Finally, applying Eq. (8) a second time gives the value for $`\alpha _s^{(4)}(\overline{\mu })`$. The described procedure is performed automatically in rhad using the subroutine rundecalpha (see App. D).

The evaluation of $`\overline{m}_Q^{(n_f)}(\mu )`$ from the input $`\overline{m}_Q(\overline{m}_Q)`$ proceeds completely analogously. The RGE that governs the running of the quark masses in the $`\overline{\text{MS}}`$ scheme reads:

$$\begin{array}{cc}\hfill \mu ^2\frac{d}{d\mu ^2}\overline{m}_Q^{(n_f)}(\mu )& =\overline{m}_Q^{(n_f)}(\mu )\gamma _m^{(n_f)}\left(\alpha _s^{(n_f)}(\mu )\right)=\overline{m}_Q^{(n_f)}(\mu )\underset{i0}{}\gamma _{m,i}^{(n_f)}\left(\frac{\alpha _s^{(n_f)}(\mu )}{\pi }\right)^{i+1}.\hfill \end{array}$$
(10)
Combining Eqs. (8) and (10) leads to

$$\begin{array}{cc}\hfill \overline{m}_Q^{(n_f)}(\mu )& =\overline{m}_Q^{(n_f)}(\mu _0)\frac{c(\alpha _s(\mu )/\pi )}{c(\alpha _s(\mu _0)/\pi )},\hfill \end{array}$$
(11)
with

$$\begin{array}{cc}\hfill c(x)& =x^{c_0}\{1+(c_1b_1c_0)x+\frac{1}{2}[(c_1b_1c_0)^2+c_2b_1c_1+b_1^2c_0b_2c_0]x^2\hfill \\ & +[\frac{1}{6}(c_1b_1c_0)^3+\frac{1}{2}(c_1b_1c_0)(c_2b_1c_1+b_1^2c_0b_2c_0)\hfill \\ & +\frac{1}{3}(c_3b_1c_2+b_1^2c_1b_2c_1b_1^3c_0+2b_1b_2c_0b_3c_0)]x^3+𝒪(x^4)\},\hfill \\ \hfill c_i& \frac{\gamma _{m,i}^{(n_f)}}{\beta _0^{(n_f)}},b_i\frac{\beta _i^{(n_f)}}{\beta _0^{(n_f)}}.\hfill \end{array}$$
(12)
If $`R(s)`$ is evaluated at $`n`$-loop order, the curly bracket of Eq. (12) has to be evaluated up to, and including, the term $`x^{n1}`$.

The matching between $`\overline{m}_Q^{(n_f)}`$ and $`\overline{m}_Q^{(n_f1)}`$ is determined through the equation

$$\begin{array}{c}\hfill \overline{m}_Q^{(n_f1)}=\zeta _m\overline{m}_Q^{(n_f)},\end{array}$$
(13)
where the function $`\zeta _m`$ can again be found in App. E.4. The procedure for the evaluation of $`\overline{m}_Q^{(n_f)}(\mu )`$ from the input $`\overline{m}_Q(\overline{m}_Q)`$ in rhad is called rundecmass (see App. D).

Note that the evaluation of $`\overline{m}_Q^{(n_f)}(\mu )`$ is only required if the user decides to use $`\overline{\text{MS}}`$ masses for the evaluation of $`R(s)`$. In this case, in parameters.f (see App. C),

```
      lmsbar = .true.

```

and the values for massc, massb, masst have to be set to the scale invariant masses $`\overline{m}_c(\overline{m}_c)`$, $`\overline{m}_b(\overline{m}_b)`$, $`\overline{m}_t(\overline{m}_t)`$, respectively. If one chooses to use pole masses instead, one defines

```
      lmsbar = .false.

```

and sets massc, massb, masst equal to the on-shell masses $`M_c`$, $`M_b`$, $`M_t`$, respectively. rhad will then simply use these values throughout the calculation.

There is one exception where we do not allow a choice between on-shell and $`\overline{\text{MS}}`$ mass: For “power-suppressed” corrections, i.e., double-bubble diagrams that contain a heavy secondary quark loop (see Fig. 3 $`(e)`$ with $`m_2\sqrt{s},m_1`$), the heavy quark mass is always inserted in the on-shell scheme. If the input for the heavy mass is provided in the $`\overline{\text{MS}}`$ scheme, the corresponding on-shell value is evaluated within rhad by the subroutine mms2mos (see App. D), using the proper conversion relations .

## 3 Notation and tree-level result

The tree level result for $`R_Q(s)`$ is part of every introductory course in quantum field theory. It reads

$$\begin{array}{c}\hfill R_Q^{(0)}(s)=n_c𝒱_Q^2r_V^{(0)}(s,m_Q),\text{with}r_V^{(0)}(s,M_Q)=\frac{v_Q}{2}(3v_Q^2),\end{array}$$
(14)
where $`n_c=3`$ is the number of colors, and $`𝒱_Q`$ is the electric charge of the quark $`Q`$ in units of the proton charge. I.e.,

$$\begin{array}{c}\hfill 𝒱_u=𝒱_c=𝒱_t=+\frac{2}{3},𝒱_d=𝒱_s=𝒱_b=\frac{1}{3}.\end{array}$$
(15)
$`v_Q`$ denotes the quark velocity,

$$\begin{array}{c}\hfill v_Q=\sqrt{1\frac{4M_Q^2}{s}},\end{array}$$
(16)
where $`\sqrt{s}`$ is the c.m.s. energy. Note that we define $`v_Q`$ in terms of the on-shell mass, $`M_Q`$ (see Sect. 2). This becomes relevant at higher orders, where we refer to the definition of Eq. (16).

## 4 One-loop result

Sample diagrams that contribute to the one-loop approximation of $`R(s)`$ are shown in Fig. 2. Their imaginary part has been evaluated a long time ago in the context of QED . In its most compact form, the result reads

$$\begin{array}{c}\hfill R_Q^{(1)}(s)=n_c𝒱_Q^2r_V^{(1)}(s,m_Q),\end{array}$$
(17)
with

$$\begin{array}{cc}\hfill r_V^{(1)}(s,M_Q)& =C_\mathrm{F}\frac{3v_Q^2}{2}[(1+v_Q^2)(\text{Li}_2\left(p^2\right)+2\text{Li}_2\left(p\right)\mathrm{ln}p\mathrm{ln}\frac{(1+v_Q)^3}{8v_Q^2})\hfill \\ & +3v_Q\mathrm{ln}\frac{1v_Q^2}{4v_Q}v_Q\mathrm{ln}v_Q+\frac{33+22v_Q^27v_Q^4}{8(3v_Q^2)}\mathrm{ln}\frac{1+v_Q}{1v_Q}+\frac{15v_Q9v_Q^3}{4(3v_Q^2)}],\hfill \end{array}$$
(18)
where $`p=(1v_Q)/(1+v_Q)`$ and $`C_\mathrm{F}=4/3`$. In the limit $`v_Q1`$, one obtains $`r_V^{(1)}=3C_\mathrm{F}/4`$ which, after taking into account the coupling and color factors, leads to the following expression for the QED corrections:

$$\begin{array}{c}\hfill \delta R_Q^{\mathrm{QED}}=n_c𝒱_Q^4\frac{3}{4}\frac{\alpha }{\pi }.\end{array}$$
(19)
The tree-level and one-loop functions and their correspondence in rhad are summarized in Tab. 1.

## 5 Two-loop result

Starting from two-loop order, a general analytic expression for $`R(s)`$ is no longer available. Typical diagrams are shown in Fig. 3. The full set can be obtained from these diagrams by attaching the two external photon lines to the quark lines at arbitrary points. Only certain classes of diagrams have been evaluated in closed form. Nevertheless, the analytic evaluation of different kinematical limits, combined with appropriate interpolation among these limits, have resulted in extremely accurate semi-analytical approximations to the full result. Thus, for all practical purposes, the full mass dependence is available at order $`\alpha _s^2`$.

Two main new features occur at the two-loop level that are absent at lower orders. First, there are contributions that are due to the non-Abelian character of QCD and have no correspondence in QED. These are diagrams with a three-gluon coupling (see Fig. 3$`(c)`$ and $`(d)`$; the four-gluon coupling occurs for the first time at three-loop order). The second new feature is that a second quark line can appear, with the effect that the result may depend on two different quark masses $`m_1`$ and $`m_2`$ (see Fig. 3$`(e)`$).

It is clear that the simplest set of diagrams to evaluate are the ones with self-energy insertions in a gluon propagator, Fig. 3$`(d)`$, $`(e)`$. In fact, the cases with massless insertions (i.e., gluons and massless quarks) are known in closed analytical form, for general values of $`s`$ and $`m_1`$ . The same is true for $`m_1=0`$, $`m_20`$  (labels in accordance with Fig. 3$`(e)`$). The case $`m_1=m_2`$ is known in terms of a two-dimensional integral representation . For $`m_1m_2`$, the coefficient of the leading term in an expansion around $`s/m_2^20`$ is known in closed form . For $`m_1=0`$ , this term was shown to approximate the exact $`s/m_2^2`$ dependence  extremely well, even up to the threshold $`s=4m_2^2`$.

The diagrams without self-energy insertions on gluon lines have not yet been evaluated in closed form for general values of $`m_Q`$ and $`s`$. The massless limit has been known for quite some time . Subsequently, mass corrections in the high energy limit have been computed: The $`m_Q^2/s`$ terms can be obtained by a simple Taylor expansion of the integrand ; the $`m_Q^4/s^2`$ terms were derived from the $`\overline{\text{MS}}`$ operator product expansion of $`\mathrm{\Pi }(q^2)`$, combined with renormalization group relations ; and the evaluation of the higher order terms (up to $`m_Q^{12}/s^6`$) was achieved by systematically expanding the Feynman integrals .

However, one should note that the convergence of such a high-energy expansion is not guaranteed below the highest threshold of the diagrams under consideration. For example, the non-planar diagram in Fig. 3$`(b)`$ has a four-quark cut at $`q^2=(4m_Q)^2`$, meaning that the expansion around $`m_Q^2/s0`$ formally converges only above $`s=(4m_Q)^2`$. Nevertheless, in practice one often observes also satisfactory convergence below the threshold .

A result that is valid in the full kinematic range was obtained in Ref.  (see also Ref. ). To this aim, the high-energy limit was combined with up to eight moments of the polarization function obtained through an expansion for $`q^20`$. This information was combined with additional input from the threshold region around $`q^2=4m_Q^2`$, using a conformal mapping and Padé approximation in order to arrive at a three-loop expression for $`\mathrm{\Pi }(q^2)`$. For all known practical applications, this approximate result is equivalent to an analytic expression for $`\mathrm{\Pi }(q^2)`$ and, after taking the imaginary part, for $`R(s)`$.

Let us now explicitely parameterize the two-loop contribution to $`R(s)`$. We can write

$$\begin{array}{cc}\hfill R_Q^{(2)}(s)=n_c𝒱_Q^2[& r_V^\mathrm{F}(s,m_Q)+r_V^\mathrm{A}(s,m_Q)+\underset{q}{}r_V^{\mathrm{db}}(s,m_Q,m_q)],\hfill \end{array}$$
(20)
where the sum runs over all quark flavors. $`r_V^\mathrm{F}C_\mathrm{F}^2`$ denotes contributions from Abelian diagrams, $`r_V^\mathrm{A}C_\mathrm{A}C_\mathrm{F}`$ comes from non-planar and non-Abelian diagrams, and $`r_V^{\mathrm{db}}C_\mathrm{F}T`$ (“double-bubble”) arises from diagrams with two separate quark lines. $`C_\mathrm{F}=4/3`$, $`C_\mathrm{A}=3`$ and $`T=1/2`$ are color factors of QCD. For later reference, it is convenient to introduce additional functions that refer to specific limits of $`r_V^{\mathrm{db}}`$:

$$\begin{array}{cc}\hfill r_V^{\mathrm{db}}(s,m_Q,0)& r_V^\mathrm{L}(s,m_Q),r_V^{\mathrm{db}}(s,m_Q,m_Q)r_V^\mathrm{T}(s,m_Q).\hfill \end{array}$$
(21)
See also Tab. 2 for more details. Note that $`r_V^{\mathrm{db}}`$ introduces contributions from secondary quarks $`q`$ other than $`Q`$, both virtual and real. They arise through the splitting of gluons into a quark–anti-quark pair.

In the massless case, one has the simple expressions

$$\begin{array}{cc}\hfill r_V^\mathrm{F}(s,m_Q=0)& =C_\mathrm{F}^2\left(\frac{3}{32}\right),\hfill \\ \hfill r_V^\mathrm{A}(s,m_Q=0)& =C_\mathrm{A}C_\mathrm{F}\left(\frac{123}{32}\frac{11}{4}\zeta _3\frac{11}{16}L_{s\mu }\right),\hfill \\ \hfill r_V^{\mathrm{db}}(s,m_1=0,m_2=0)& =C_\mathrm{F}T\left(\frac{11}{8}+\zeta _3+\frac{1}{4}L_{s\mu }\right).\hfill \end{array}$$
(22)
where $`L_{s\mu }\mathrm{ln}(s/\mu ^2)`$ and $`\zeta _31.20206`$. The expressions in the massive case are given in App. E.1.

In the following, the expressions adopted for the individual contributions are described, and their origin is given.

* For $`r_V^\mathrm{F}(s,m_Q)`$ and $`r_V^\mathrm{A}(s,m_Q)`$ we use the full mass dependence as given by the Padé results of Eqs. (36), (39.
* For the double-bubble contribution $`r_V^{\mathrm{db}}(s,m_1,m_2)`$, we have to distinguish several cases:

+ $`m_1=m_2=m_Q`$: For $`\sqrt{s}4m_Q`$, we use the analytic formula as given in Eq. (44. For $`\sqrt{s}>4m_Q`$, the full mass dependence is known in terms of a two-fold integral representation, see Eq. (46. For the sake of speed, however, we use the high energy expansion of Eq. (48. Numerical differences to the integral representation are completely negligible. Nevertheless, the user may switch to the integral formula by setting lrvctexp = .false. in rhad.
    + $`m_2<m_1`$: We apply an expansion for $`m_2^2m_1^2(<s/4)`$:

    1. At $`𝒪(m_2^0)`$, the full $`m_1`$-dependence is used in the form of the analytic result of Eq. (41) (for $`x=\mathrm{L}`$.
        2. At $`𝒪(m_2^2)`$, we neglect the $`m_1`$ dependence. It turns out that this contribution vanishes which can be seen in Eq. (52).
        3. Higher orders in $`m_2`$ are neglected.
    + $`(2m_1)^2<s<(2m_2)^2`$: If $`m_1=0`$, we use the analytical result of Eq. (49, including the full dependence on $`m_2`$. If $`m_10`$, we use the leading term in $`1/m_2^2`$, keeping its full $`m_1^2`$ dependence, see Eq. (54). This expression has been obtained in Ref. .
    + $`(2m_1)^2<(2m_2)^2<s`$: If $`m_1=0`$, we use the analytical result of Eq. (49, including the full dependence on $`m_2`$. If $`m_10`$, we apply an expansion in the limit $`m_1^2m_2^2(<s/4)`$:

    1. At $`𝒪(m_1^0)`$, we use the analytical result of Eq. (49 for the full $`m_2`$ dependence.
        2. At $`𝒪(m_1^2)`$, we neglect all effects from non-zero $`m_2`$. The corresponding expression for $`r_V^{\mathrm{db}}(s,m_1,0)`$ is given in Eq. (LABEL:eq::rvdbm1).
        3. Higher orders in $`m_1`$ are neglected.

Note that we set the ratios $`m_c^2/m_t^2`$ and $`m_b^2/m_t^2`$ to zero and consider mass corrections due to the presence of an additional light quark only for $`c`$-quark effects in $`b`$-quark production.

For the sake of completeness, let us remark that also the two-loop corrections of order $`\alpha \alpha _s`$ and order $`\alpha ^2`$ for massless quarks  are known. Numerically, these contributions are very small. Nevertheless, we include the mixed corrections of order $`\alpha \alpha _s`$ into rhad, which modifies Eq. (19) to

$$\begin{array}{c}\hfill \delta R_Q^{\mathrm{QED}}=n_c𝒱_Q^4\frac{3}{4}\frac{\alpha }{\pi }\left(1\frac{1}{3}\frac{\alpha _s(\mu )}{\pi }\right).\end{array}$$
(23)
Higher order QED or mixed QED/QCD effects are numerically irrelevant and will be neglected .

## 6 Three-loop result

At order $`\alpha _s^3`$, a new class of diagrams contributes to $`R_{\mathrm{had}}`$, the so-called “singlet diagrams.” They differ from the non-singlet diagrams in the sense that, in the singlet case, the external currents are not connected by a common quark line. A typical diagram is shown in Fig. 4$`(d)`$. The non-singlet diagrams (see Fig. 4$`(a)`$$`(c)`$) are numerically dominant, but rhad includes both contributions. We will describe them in more detail in what follows.

### 6.1 Non-singlet contributions

The knowledge of $`R(s)`$ at three-loop level is restricted to the high energy limit. The massless limit was obtained for the first time in Ref.  and later confirmed through an independent calculation in Ref. .

Mass corrections are known up to $`𝒪(m_Q^4/s^2)`$: The quadratic terms were obtained from renormalization group identities , while the quartic mass terms were evaluated by applying the methods of Ref.  at three-loop level .

Let us briefly discuss how we treat mass effects from massive inner quark loops. At three-loop level, there can be three different quark types in one diagram, with, say $`m_1>m_2>m_3`$. The approximation that we apply is to take all terms including $`m_1^4`$ and $`m_2^2`$, neglecting mixed terms of order $`m_1^2m_2^2`$. Mass effects from $`m_3`$ are also neglected. In addition, we do not consider power-suppressed terms, i.e. diagrams that contain a quark $`q`$ with $`s<s_q^{\mathrm{thr}}`$.

With these approximations in mind, we can write the three-loop contribution as

$$\begin{array}{cc}\hfill R_Q^{(3)}(s)& =n_c𝒱_Q^2\left[r_V^{(3)}(s,m_Q,n_f)+\underset{qQ}{}\delta r_0^{(3)}(s,m_q,n_f)\right].\hfill \end{array}$$
(24)
$`\delta r_0^{(3)}(s,m_q,n_f)`$ denotes the mass corrections due to massive inner quark loops different from $`Q`$. According to the approximations above, it vanishes for $`s<s_q^{\mathrm{thr}}`$ and also for $`m_q=0`$.

Note that $`n_f`$ itself depends on the c.m.s. energy. Assume, for example, that $`s<s_c^{\mathrm{thr}}`$. Then we have $`n_f=3`$ and

$$\begin{array}{cc}\hfill R_Q^{(3)}(s)& =n_c𝒱_Q^2r_V^{(3)}(s,0,3),Q=u,d,s.\hfill \end{array}$$
(25)
For $`s_c^{\mathrm{thr}}<s<s_b^{\mathrm{thr}}`$, we have $`n_f=4`$ and thus the non-zero contributions are

$$\begin{array}{cc}\hfill R_Q^{(3)}(s)& =n_c𝒱_Q^2\left[r_V^{(3)}(s,0,4)+\delta r_0^{(3)}(s,m_c,4)\right],Q=u,d,s,\hfill \\ \hfill \text{and}R_c^{(3)}(s)& =n_c𝒱_c^2r_V^{(3)}(s,m_c,4).\hfill \end{array}$$
(26)
For $`s_b^{\mathrm{thr}}<s<s_t^{\mathrm{thr}}`$, it is $`n_f=5`$ and

$$\begin{array}{cc}\hfill R_Q^{(3)}(s)& =n_c𝒱_Q^2\left[r_V^{(3)}(s,0,5)+\delta r_0^{(3)}(s,m_c,5)+\delta r_0^{(3)}(s,m_b,5)\right],Q=u,d,s,\hfill \\ \hfill R_c^{(3)}(s)& =n_c𝒱_c^2\left[r_V^{(3)}(s,m_c,5)+\delta r_0^{(3)}(s,m_b,5)\right],\hfill \\ \hfill R_b^{(3)}(s)& =n_c𝒱_b^2\left[r_V^{(3)}(s,m_b,5)+\delta r_0^{(3)}(s,m_c,5)\right].\hfill \end{array}$$
(27)
Finally, for $`s>s_t^{\mathrm{thr}}`$, all masses except for $`m_t`$ are neglected, and, with $`n_f=6`$, we have

$$\begin{array}{cc}\hfill R_Q^{(3)}(s)& =n_c𝒱_Q^2\left[r_V^{(3)}(s,0,6)+\delta r_0^{(3)}(s,m_t,6)\right],Q=u,d,c,s,b,\hfill \\ \hfill R_t^{(3)}(s)& =n_c𝒱_t^2r_V^{(3)}(s,m_t,6).\hfill \end{array}$$
(28)
As mentioned before, in Eq. (27) it is understood that all $`c`$-quark mass effects are included only at $`𝒪(m_c^2)`$.

The expressions for $`r_V^{(3)}(s,m_Q,n_f)`$ and $`\delta r_0^{(3)}(s,m_Q,n_f)`$ are listed in App. E.2.

### 6.2 Singlet contributions

As already mentioned above, singlet diagrams are characterized by the fact that each of the external photons couples to a separate quark line (see, e.g., Fig. 4$`(d)`$). The singlet contribution arises first at three-loop level; at two-loop level, where the two quark lines are connected by only two gluons, it is zero due to Furry’s theorem.

We write the singlet contribution to $`R(s)`$ in the following way:

$$\begin{array}{cc}\hfill R_{\mathrm{sing}}(s)& =n_c\underset{q,q^{}}{}𝒱_q𝒱_q^{}r_{V,\mathrm{sing}}^{(3)}(s,m_q,m_q^{}),\hfill \end{array}$$
(29)
where the sum over $`q`$ and $`q^{}`$ runs over all active quark flavors. E.g., for $`\sqrt{s}=12`$ GeV it includes all quarks except for the top quark. Again, only the high energy expansion is known. Quadratic mass terms are absent, and the expression including quartic mass terms reads:

$$\begin{array}{cc}\hfill r_{V,\mathrm{sing}}^{(3)}(s,m_1,m_2)& =\frac{55}{216}\frac{5}{9}\zeta _3+\frac{m_1^4}{s^2}\left(\frac{10}{9}+\frac{25}{3}\zeta _3\right)+\frac{m_2^4}{s^2}\left(\frac{10}{9}+\frac{25}{3}\zeta _3\right).\hfill \end{array}$$
(30)
Note that, if $`m_1m_2`$, we neglect the lighter of the two masses. Since this is the lowest order where singlet terms occur, Eq. (30) is the same in the $`\overline{\text{MS}}`$ and in the on-shell mass scheme.

## 7 Four-loop result

The knowledge of the four-loop contribution to $`R(s)`$ is very limited. Only contributions of diagrams with two and three quark-loop insertions have been calculated. Sample diagrams are shown in Fig. 5$`(c)`$ and $`(d)`$. The mass effects are evaluated up to order $`m_Q^2/s`$, so that the results are strictly only valid in the high-energy limit.

We write the four-loop contribution as

$$\begin{array}{cc}\hfill \overline{R}_Q^{(4)}(s)& =n_c𝒱_Q^2\overline{r}_V^{(4)}(s,m_Q),\hfill \end{array}$$
(31)
with

$$\begin{array}{cc}\hfill \overline{r}_V^{(4)}(s,m_Q)& =\overline{r}_0^{(4),0\mathrm{L}}(s)+n_f\overline{r}_0^{(4),1\mathrm{L}}(s)+n_f^2r_0^{(4),2\mathrm{L}}(s)+n_f^3r_0^{(4),3\mathrm{L}}(s)\hfill \\ & +\frac{m_Q^2}{s}\left(\overline{r}_{V,2}^{(4),0\mathrm{L}}(s)+n_f\overline{r}_{V,2}^{(4),1\mathrm{L}}(s)+n_f^2r_{V,2}^{(4),2\mathrm{L}}(s)+n_f^3r_{V,2}^{(4),3\mathrm{L}}(s)\right),\hfill \end{array}$$
(32)
where the bar indicates that the expressions are estimates rather than approximations or even exact results. $`n_f`$ is the number of active flavors, varying with the c.m.s. energy as discussed in Sect. 2.

$`r_0^{(4),3\mathrm{L}}`$ is the renormalon-type contribution with three massless quark insertions on the gluon propagator, see Fig. 5$`(d)`$. It has been evaluated in Ref. , where the general structure of the terms of order $`\alpha _s(\alpha _sn_f)^n`$ was derived. Recently, the analytic expression for $`r_0^{(4),2\mathrm{L}}`$ became available . It required the evaluation of massless four-loop two-point functions in combination with the method of Ref.  (see also Ref. ) to derive the imaginary part of the five-loop contributions as shown in Fig. 5$`(c)`$. The same method has been used in combination with the technique derived in Ref.  for the computation of the quadratic mass corrections $`r_{V,2}^{(4),2\mathrm{L}}`$ and $`r_{V,2}^{(4),3\mathrm{L}}`$ .

Estimates for the full massless four-loop result have been known before and are still a subject of interest . We determine $`\overline{r}_0^{(4)}(s)`$ and $`\overline{r}_0^{(4),\mathrm{L}}(s)`$ by subtracting the known $`n_f^2`$ and $`n_f^3`$ contributions from the estimates of Ref.  at $`n_f=1,\mathrm{},6`$, and fitting the resulting six “data points” by a linear function in $`n_f`$. One finds, for $`\mu ^2=s`$,

$$\begin{array}{cc}\hfill \overline{r}_0^{(4),0\mathrm{L}}(s)& =1.8610^2,\overline{r}_0^{(4),1\mathrm{L}}(s)=21.3,\hfill \\ \hfill r_0^{(4),2\mathrm{L}}(s)& =7.9710^1,r_0^{(4),3\mathrm{L}}(s)=2.1510^2.\hfill \end{array}$$
(33)
The logarithmic contributions follow from the lower order terms through renormalization group invariance and are collected in App. E.3. Once the exact results for $`r_0^{(4),0\mathrm{L}}`$ and $`r_0^{(4),1\mathrm{L}}`$ become available, the approximate ones can easily be replaced in rhad.

Using the same method for the quadratic mass terms in combination with the estimates of Ref.  we obtain, in the $`\overline{\text{MS}}`$ scheme

$$\begin{array}{cc}\hfill \overline{r}_{V,2}^{(4),0\mathrm{L}}(s)& =7.1110^3,\overline{r}_{V,2}^{(4),1\mathrm{L}}(s)=1.4310^3,\hfill \\ \hfill r_{V,2}^{(4),2\mathrm{L}}(s)& =49.1,r_{V,2}^{(4),3\mathrm{L}}(s)=0.204.\hfill \end{array}$$
(34)
In App. E.3, the corresponding result for on-shell quark masses is listed together with the logarithmic contributions.

For completeness, the function implemented in rhad is listed in Tab. 4. Note that at $`𝒪(\alpha _s^4)`$, the corrections analoguous to $`\delta r_0^{(3)}`$ are neglected.

## 8 Evaluating $`𝑹\mathbf{(}𝒔\mathbf{)}`$

In this section we describe how rhad evaluates $`R(s)`$ in the individual energy regions. Let us first recall that $`s_Q^{\mathrm{thr}}`$ defines the lowest value of the c.m.s. energy squared, $`s`$, at which the perturbative treatment of $`Q\overline{Q}`$ production is allowed; $`R_Q(s)`$ is defined to be zero for $`s<s_Q^{\mathrm{thr}}`$. The physical threshold for $`Q\overline{Q}`$ production is at $`s_Q^{\mathrm{low}}`$; the user is advised, however, to disregard the results of rhad in the region between $`s_Q^{\mathrm{low}}`$ and $`s_Q^{\mathrm{thr}}`$ ($`Q=c,b,t`$).

The evaluation of the number of active flavors $`n_f`$, the strong coupling constant $`\alpha _s^{(n_f)}(\mu )`$, and, if required, the $`\overline{\text{MS}}`$ quark masses $`\overline{m}_Q^{(n_f)}(\mu )`$, from the input values $`s`$, $`\mu `$, $`\alpha _s^{(5)}(M_Z)`$, and $`\overline{m}_Q(\overline{m}_Q)`$, has been described in Sect. 2. Let us now look in more detail at the specific contributions that enter $`R(s)`$ at certain values of $`s`$.

For $`s<s_c^{\mathrm{low}}`$, i.e., below the production threshold for two $`D_0`$ mesons, we have $`n_f=3`$. The charm, bottom and top quark masses are decoupled and thus their contribution goes like $`\alpha _s^2s/m_q^2`$ ($`q=c,b,t`$), see Eqs. (49,54). All other correction terms are evaluated in the massless limit.

In the region $`s_c^{\mathrm{thr}}<s<s_b^{\mathrm{low}}`$ one has $`n_f=4`$. $`\alpha _s^{(4)}(\mu )`$ is evaluated from the input $`\alpha _s^{(5)}(M_Z)`$, and, in the $`\overline{\text{MS}}`$ mass scheme, $`\overline{m}_c^{(4)}(\mu )`$ is obtained from the input quantity $`\overline{m}_c(\overline{m}_c)`$ (see Sect. 2). In the on-shell mass scheme, the input value $`M_c`$ is used unchanged throughout the calculation. The full charm quark mass dependence (inasmuch it is known) is taken into account. Bottom and top quark are decoupled and enter through the power-suppressed terms, Eqs. (49,54), at order $`\alpha _s^2`$.

At c.m.s. energies above $`s_b^{\mathrm{thr}}`$ and below $`s_t^{\mathrm{low}}`$, the situation is similar to the previously described region, but with $`n_f=5`$. Note, however, that in the $`\overline{\text{MS}}`$ scheme not only $`\overline{m}_b^{(5)}(\mu )`$ is needed for the numerical evaluation, but also $`\overline{m}_c^{(5)}(\mu )`$. It is obtained from the input $`\overline{m}_c(\overline{m}_c)`$ using the routine rundecmass (see Sect. D).

For $`s>s_t^{\mathrm{thr}}`$, $`\alpha _s^{(6)}(\mu )`$ and $`\overline{m}_t^{(6)}(\mu )`$ are needed for the evaluation of $`R(s)`$. At these energies, no mass corrections from lighter quarks are considered.

Tab. 5 and 6 show the results for $`R(s)`$ with the default settings (see Sect. C), for different values of the c.m.s. energy. The individual rows correspond the different orders of perturbation theory. In the case of $`\sqrt{s}=5`$ GeV and $`\sqrt{s}=12`$ GeV, we also show the partial contributions to $`R(s)`$ arising from massive quarks. The second column contains the value of $`\alpha _s^{(n_f)}(\sqrt{s})`$, with the appropriate $`n_f`$.

Tab. 7 contains the results for $`\alpha _s(s)`$, $`R(s)`$, $`R_c(s)`$ and $`R_b(s)`$ for the energy region between $`\sqrt{s}=1.8`$ GeV and $`\sqrt{s}=20`$ GeV, sparing out those parts where perturbation theory is not applicable. We adopted the default values of rhad given in App. C with iord = 4. In Fig. 6 the data for $`R(s)`$ (full curve) are shown in graphical form. The shaded bands mark the uncertainty of $`R(s)`$ induced by the variation of $`\alpha _s(M_Z)`$, $`\mu `$, $`M_c`$ and $`M_b`$ in the range

$`\alpha _s(M_Z)`$ $`=`$ $`0.118\pm 0.003,`$
$`M_c`$ $`=`$ $`(1.65\pm 0.15)\text{GeV},`$
$`M_b`$ $`=`$ $`(4.75\pm 0.2)\text{GeV},`$
$`\mu `$ $``$ $`[{\displaystyle \frac{\sqrt{s}}{2}},2\sqrt{s}].`$ (35)
## 9 A typical program

It is instructive to look at a typical program that evaluates $`R(s)`$. We set the c.m.s. energy to $`\sqrt{s}=12`$ GeV, and split the result according to the contributions from the individual quarks. The occuring functions will be described in more detail in App. B, but their functionality should be clear from the following program.

#### Input.

The corresponding fortran program would look as follows.

```
      program example

      implicit real*8(a-h,m-z)
      implicit integer(i,j)
      implicit character*60(k)
      implicit logical(l)

      include ’common.f’

      sqrts = 12.d0
      scms  = sqrts*sqrts

      call parameters(scms)
      call init(scms)

      rall = rhad(scms)
      ru = ruqrk(scms)
      rd = rdqrk(scms)
      rs = rsqrk(scms)
      rc = rcqrk(scms)
      rb = rbqrk(scms)
      rt = rtqrk(scms)
      rsg = rsinglet(scms)
      rem = rqed(scms)

      print*,’R_had  = ’,rall
      print*,’R_u    = ’,ru
      print*,’R_d    = ’,rd
      print*,’R_s    = ’,rs
      print*,’R_c    = ’,rc
      print*,’R_b    = ’,rb
      print*,’R_t    = ’,rt
      print*,’R_sing = ’,rsg
      print*,’R_QED  = ’,rem

      end

```

#### Output.

If lverbose is set to .true. the output reads:

```

            rhad.f -- Version 1.00
 by Robert Harlander and Matthias Steinhauser
                December 2002

Order of calculation:   4

Scales:
  sqrt(s)  =    12.000 GeV
  mu       =    12.000 GeV
  thrc     =     4.800 GeV
  thrb     =    11.200 GeV
  thrt     =   360.000 GeV

Number of active flavors:
  nf  =  5

Coupling constants:
  alpha_QED    =  1 / ( 137.036 )
  alpha_s(Mz)  =  0.1180             [ Mz =  91.1876 GeV ]
  alpha_s(mu)  =  0.1709360043

Quark masses:
  M_c       =    1.65 GeV
  M_b       =    4.75 GeV
  M_t       =  175.00 GeV

General switches (F=False, T=True):
  only massless terms    : F
  power suppressed terms : T
  QED corrections        : T
  singlet contributions  : T
  alphas^3 m^2 included  : T
  alphas^3 m^4 included  : T
  alphas^4 m^2 included  : T
 -- end of parameters --
 R_had  =   3.9049832
 R_u    =   1.40765446
 R_d    =   0.351913615
 R_s    =   0.351913615
 R_c    =   1.41577367
 R_b    =   0.375554797
 R_t    =   0.
 R_sing =  -4.42977631E-05
 R_QED  =   0.00221733842

```

If lverbose = .false., only the lines after “-- end of parameters --” are printed. The parameter list displays the most important initializations, like the value of $`\alpha _s^{(5)}(M_Z)`$ and the renormalization scale $`\mu `$. It also contains values for parameters that were derived from these settings, like the number of active flavors $`n_f`$, and the strong coupling constant $`\alpha _s^{(n_f)}(\mu )`$.

## 10 Conclusions

In this paper we have discussed the perturbative corrections to the inclusive cross section $`\sigma (e^+e^{}\text{hadrons})`$ and collected the analytical and semi-analytical formulas. This includes the full mass dependence up to order $`\alpha _s^2`$, the expansion up to quartic mass terms at order $`\alpha _s^3`$ and the quadratic mass corrections at order $`\alpha _s^4`$. Furthermore, the running and decoupling formalism necessary for a consistent evaluation of the strong coupling constant and the $`\overline{\text{MS}}`$ quark masses has been presented.

The main subject of this paper, however, is a description of the fortran program rhad which allows for the evaluation of $`R(s)`$, including all currently available radiative corrections. The program is straightforward to use, in particular if the default parameter set is adopted. Variations of the physical input values should be done such that they still resemble the physical case. The modularity of rhad allows for a simple extension once new corrections to the theoretical predictions for $`R(s)`$ become available. In conclusion, rhad can easily be used to compute, for example, the perturbative parts of the hadronic contributions to the running of the electromagnetic coupling, and the anomalous magnetic moment of the muon.

Acknowledgments

First of all, we are indebted to Hans Kühn for suggesting this project and for his valuable advice at all stages of its development. It is a pleasure to thank Thomas Teubner for his help and patience in debugging rhad. The comparison to his private code was an essential ingredient for the completion of this work. Furthermore, we thank him and Marek Jeżabek for providing fortran routines which have been used in preliminary versions of rhad. Thorsten Seidensticker kindly allowed us to include the unpublished analytical results of $`r_V^{\mathrm{db}}(s,m_1,m_2)`$ in the limit $`(2m_1)^2<s<(2m_2)^2`$. Furthermore, we thank the authors of Refs.  for making the manuscripts accessible to us before publication.

Appendix

## Appendix A Installation

The distribution of rhad contains the following files:

```
Examples  example.f  makefile      r012.f  rhad.f   vegas-rhad.f
common.f  funcs.f    parameters.f  r34.f   runal.f

```

Examples is a directory that contains programs to reproduce Tables 5, 6, and 7 of this paper. One example program is kept in the main directory, example.f. Its listing in shown in Sect. 9. It can be compiled by simply calling GNU make:

```
> gmake prog=example

```

The executable is named xexample.

If GNU make is not available, the fortran files can be compiled individually using

```
f77 -c -o <file>.o <file>.f

```

and then linked using

```
f77 -o xexample *.o

```

## Appendix B Basic functions of rhad

Unless indicated otherwise, we use the following implicit type specifications in rhad:

```
      IMPLICIT  REAL*8(A-H,M-Z)
      IMPLICIT  INTEGER(I,J)
      IMPLICIT  CHARACTER*60(K)
      IMPLICIT  LOGICAL(L)

```

The functions that encode the different loop orders to $`R(s)`$ have already been discussed in Sects. 37, and in particular in Tables 14.

The proper sums of these contributions are taken by the following functions corresponding to $`R_Q(s)`$, $`R_{\mathrm{sing}}(s)`$ and $`\delta R_Q^{\mathrm{QED}}(s)`$ (photon exchange only):

$`R_d(s)`$, $`R_u(s)`$, $`R_s(s)`$, $`R_c(s)`$, $`R_b(s)`$, $`R_t(s)`$

$`R_{\mathrm{sing}}(s)`$

$`_Q\delta R_Q^{\mathrm{QED}}(s)`$

The full $`R`$-ratio is obtained by calling rhad(s), which adds the contributions from the individual quarks, depending on the c.m.s. energy.

The general structure of a fortran program from which the above functions are called is as follows:

```
      program <name of program>
      <declaration of variables>
      <definition of s, e.g.:>
      include ’common.f’
      scms = 11.5d0**2
      call parameters(scms)
      call init(scms)
      <call of function to compute R(s)>
      end

```

It is important to call the subroutines parameters and init (in this order), which define the parameters and evaluate $`\alpha _s(\mu )`$ and, if required, evolve and convert the masses. An explicit example program has been given in Sect. 9.

## Appendix C The subroutine parameters

### C.1 Description of the parameters

The file parameters.f contains the subroutine parameters and collects the most important variables that can be adjusted by the user. parameters.f should be understood as the input file for rhad. In Appendix C.2, our default settings are listed. Let us stress that, although many of the implemented formulas are valid over a large parameter range, one should vary the input parameters only within a sensible region about their default values. Unreasonable settings, like inverted mass hierarchies or the like, may lead to inconsistent results.

The following list describes the variables defined in parameters.f:

Print values of parameters (see Sect. 8)

Output unit for parameter list.

Order of the calculation. Values: iord=0 (Born) to iord=4. It also governs the order at which the running coupling and the $`\overline{\text{MS}}`$ masses are evaluated.

$`\alpha _s^{(5)}(M_Z)`$. Input value for the strong coupling constant (see Sect. 2).

If .true., evaluate $`R(s)`$ by using $`\overline{\text{MS}}`$ masses (see Sect. 2).

Initial values for $`m_c`$, $`m_b`$, $`m_t`$. If lmsbar = .true.: scale invariant mass $`\overline{m}_q(\overline{m}_q)`$; otherwise: on-shell mass $`M_q`$ (see Sect. 2).

Renormalization scale $`\mu `$ at which $`\alpha _s`$, the $`\overline{\text{MS}}`$ quark masses and $`R(s)`$ are evaluated (see Sect. 2).

$`\mu _c`$, $`\mu _b`$, $`\mu _t`$. Matching scales for the decoupling of the heavy quarks (see Sect. 2).

If .true., include the lowest order QED and, if iord.ge.2, the mixed QED/QCD corrections (see Eqs. (19) and (23)).

If .true., neglect all quark masses.

If .true., include power-suppressed terms (see Eqs. (51),(54)).

If .true., include $`\alpha _s^3m_Q^2`$ terms (see Sect. E.2).

If .true., include $`\alpha _s^3m_Q^4`$ terms (see Sect. E.2).

If .true., include singlet terms (see Sect. 6.2).

If .true., include $`\alpha _s^4m_Q^2`$ terms (see Sect. 7).

$`s_c^{\mathrm{thr}}`$, $`s_b^{\mathrm{thr}}`$, $`s_t^{\mathrm{thr}}`$ (see Sect. 1 and 8).
Thresholds for open (perturbative) $`c`$-, $`b`$-, and $`t`$-quark pair production.

$`s_c^{\mathrm{low}}`$, $`s_b^{\mathrm{low}}`$, $`s_t^{\mathrm{low}}`$ (see Sect. 1 and 8).
Perturbation theory is not applicable for $`s_Q^{\mathrm{low}}ss_Q^{\mathrm{thr}}`$, $`Q\{c,b,t\}`$. Note that these three parameters have no influence on the results; their only task is to remind the user of the problematic energy regions in the perturbative approach.

minimal allowed value for the c.m.s. energy (see Sect. 1).

Note that not all parameters are independent: for example, if lmassless = .true., la3m2, la3m4, and la4m2 will automatically be set to .false..

### C.2 Default settings

The default settings in parameters.f are as follows:

```
      subroutine parameters(scms)
c..
c..   User-defined parameters.
c..
      implicit  real*8(a-h,m-z)
      implicit  integer(i,j)
      implicit  character*60(k)
      implicit  logical(l)

      include ’common.f’

c..   verbose mode:
      lverbose = .true.

c..   output unit for parameter list (6 = STDOUT)
      iunit = 6

c..   order or calculation:
      iord  = 4

c..   strong coupling constant at scale mz (5 active flavors):
      alphasmz  = 0.118d0

c..   use MS-bar or pole quark mass?  (.true. == MS-bar mass)
      lmsbar    = .false.

c..   masses
      massc      = 1.65d0        ! charm
      massb      = 4.75d0        ! bottom
      masst      = 175.d0        ! top

c..   renormalization scale:
      mu        = dsqrt(scms)

c..   decoupling scales:
      muc       = 2.d0*massc    ! charm
      mub       = massb         ! bottom
      mut       = masst         ! top

c..   threshold for open quark production
      thrc = 4.8d0
      thrb = 11.2d0
      thrt = 2*masst+10.d0

c..   lower bound of quark threshold region
      thrclow = 3.73d0
      thrblow = 10.52d0
      thrtlow = 2*masst-10.d0

c..   minimum allowed cms energy:
      sqmin     = 1.8d0

c..   some switches
      lmassless = .false.        ! use only massless approximation
      lqed      = .true.         ! QED corrections (.true. = ON)
      lpsup     = .true.         ! include power suppressed terms
      la3m2     = .true.         ! \alpha_s^3 * m^2 terms
      la3m4     = .true.         ! \alpha_s^3 * m^4 terms
      la3sing   = .true.         ! singlet contribution at order \alpha_s^3
      la4m2     = .true.         ! \alpha_s^4 * m^2 terms

      return
      end

```

## Appendix D Implementation of running and decoupling

#### Determination of $`𝜶_𝒔\mathbf{.}`$

rhad contains several routines which cover the consistent running and decoupling of the strong coupling $`\alpha _s^{(n_f)}(\mu )`$ (in this context see also Refs. ).

$`\alpha _s^{(n_f)}(\mu _0)/\pi `$

$`\mu _0`$

$`\mu `$

$`n_f`$

number of loops used for QCD $`\beta `$ function

0=quiet, 1=verbose

$`\alpha _s^{(n_f)}(\mu )/\pi `$

Evaluates $`\alpha _s^{(n_f)}(\mu )`$ from $`\alpha _s^{(n_f)}(\mu _0)`$ by solving numerically the renormalization group equation

$`n_i`$

$`M_h`$, heavy quark mass to be decoupled

$`\mu ^{\mathrm{th}}`$

$`k_{\mathrm{ord}}`$

$`\delta \{1,1\}`$

$`\alpha _s^{(n_i)}`$

$`\alpha _s^{(n_i+\delta )}`$

Evaluates $`\alpha _s^{(n_i+\delta )}(\mu ^{\mathrm{th}})`$ from $`\alpha _s^{(n_i)}(\mu ^{\mathrm{th}})`$. The order of the matching relations is determined by $`k_{\mathrm{ord}}`$ ($`0k_{\mathrm{ord}}3`$). The matching coefficients require the pole quark mass as input.

$`\mu _0`$

$`\mu _1`$

$`\alpha _s^{(n_i)}(\mu _0)`$, where $`n_i`$ is defined in the subroutine parameters

$`\alpha _s^{(n_f)}(\mu _1)`$, where $`n_f`$ is determined with the help of the variables thrc, thrb and thrt

Determines $`\alpha _s^{(n_f)}(\mu _1)`$ from $`\alpha _s^{(n_i)}(\mu _0)`$, including decoupling and matching, by calling decalpha and runalpha. The value for $`n_i=`$infini is set in the subroutines parameters. $`n_f=`$inffin is determined according to the value of mu1 and becomes part of the common-block.

In addition there is the subroutine decalphams which works analoguously to decalpha, except that the matching coefficients are parameterized in terms of $`\overline{\text{MS}}`$ quark masses.

#### Determination of $`\overline{𝒎}_𝑸\mathbf{(}\sqrt{𝒔}\mathbf{)}`$.

$`\overline{m}_Q(\mu _0)`$

$`\alpha _s(\mu _0)/\pi `$

$`\alpha _s(\mu _f)/\pi `$

$`n_f`$

number of loops

$`\overline{m}_Q(\mu _f)`$

Evaluates $`\overline{m}_Q^{(n_f)}(\mu )`$ from $`\overline{m}_Q^{(n_f)}(\mu _0)`$ using Eq. (12).

$`\overline{m}_Q^{(n_i)}(\mu ^{\mathrm{th}})`$

$`\alpha _s^{(n_i)}(\mu ^{\mathrm{th}})`$

$`\overline{m}_h`$, $`\overline{\text{MS}}`$ value of heavy quark mass to be decoupled

$`\mu ^{\mathrm{th}}`$ decoupling scale

$`n_i`$

$`\delta \{1,1\}`$

$`k_{\mathrm{ord}}`$

$`\overline{m}_Q^{(n_i+\delta )}(\mu ^{\mathrm{th}})`$

Evaluates $`\overline{m}_Q^{(n_i+\delta )}(\mu ^{\mathrm{th}})`$ from $`\overline{m}_Q^{(n_i)}(\mu ^{\mathrm{th}})`$. The order of the matching relations is determined by $`k_{\mathrm{ord}}`$ ($`0k_{\mathrm{ord}}3`$). The matching coefficients require the $`\overline{\text{MS}}`$ quark mass as input.

$`\overline{m}_Q^{(n_i)}(\mu _0)`$

$`n_i`$

$`\mu _0`$

$`\mu _1`$

$`\overline{m}_Q^{(n_f)}(\mu _1)`$

Determines the $`\overline{\text{MS}}`$ mass $`\overline{m}_Q^{(n_f)}(\mu _1)`$ from $`\overline{m}_Q^{(n_i)}(\mu _0)`$, including decoupling and matching, by calling decmassms and runmass. $`n_f`$ is determined according to the value of mu1.

Note that in decmassms and rundecmass only those cases are implemented which are needed for rhad.

#### Pole mass versus $`\overline{\text{MS}}`$ mass.

In case the $`\overline{\text{MS}}`$ quark masses are used as input the numerical values for the pole quark masses, which are needed for certain (double-bubble) contributions as described in Sect. 2, are computed with the help of the routine mms2mos.

$`\overline{m}_Q(\overline{m}_Q)`$

number of active flavors

number of loops used for the conversion

$`M_Q`$

Determines the on-shell mass $`M_Q`$ from the scale-invariant mass $`\overline{m}_Q(\overline{m}_Q)`$, where the inloop-loop relation is used. In practice we use inloop=iord. For the conversion also $`\alpha _s(\overline{m}_Q)`$ is needed which we determine with the help of the subroutines runalpha and decalphams.

## Appendix E Analytical results

All results presented in this section are parametrized in terms of pole quark masses denoted by $`M_Q`$, $`M_1`$ or $`M_2`$. The conversion to the $`\overline{\text{MS}}`$ formulas can be performed with the help of Eq. (67).

### E.1 Analytical results for $`R_Q^{(2)}(s)`$

In this appendix we list the results for the two-loop contributions to $`R(s)`$ as they are implemented in rhad. They are either known analytically, or in terms of semi-analytical Padé approximants.

The Abelian two-loop contribution reads

$`r_V^\mathrm{F}(s,M_Q)`$ $`=`$ $`r_V^{\mathrm{F},pa}4C_\mathrm{F}r_V^{(1)}(s,M_Q),`$ (36)
where $`r_V^{\mathrm{F},pa}`$ is known in the form of Padé approximants . We quote the result based on a $`[5/5]`$-Padé approximant:

$`r_V^{\mathrm{F},pa}`$ $`=`$ $`C_\mathrm{F}^2{\displaystyle \frac{3}{4\pi }}\text{Im}\{{\displaystyle \frac{i\pi }{48z}}[(687+186z)\sqrt{11/z}+216\mathrm{ln}{\displaystyle \frac{1\sqrt{11/z}}{1+\sqrt{11/z}}}]+{\displaystyle \frac{(1+\stackrel{~}{\omega })^2}{1\stackrel{~}{\omega }}}`$ (37)
$`\times [2659.1447467995658619093787+994.90626341783737509012835\stackrel{~}{\omega }`$
$`2742.9905230464882915790412\stackrel{~}{\omega }^2779.40067776480109262194206\stackrel{~}{\omega }^3`$
$`+562.98408428934139747541229\stackrel{~}{\omega }^4+132.845153672444424204018458\stackrel{~}{\omega }^5]`$
$`\times [121.9768095064564342458129488+11.5590492918470301261688922\stackrel{~}{\omega }`$
$`119.511603603837393164213484\stackrel{~}{\omega }^21.04113850838776740702415214\stackrel{~}{\omega }^3`$
$`+18.4653468020861357442841578\stackrel{~}{\omega }^4+\stackrel{~}{\omega }^5]^1\},`$
and

$`z`$ $`=`$ $`{\displaystyle \frac{1}{1v_Q^2}},v_Q=\sqrt{1{\displaystyle \frac{4M_Q^2}{s}}},\stackrel{~}{\omega }={\displaystyle \frac{2}{z}}1+2i\sqrt{{\displaystyle \frac{1}{z}}{\displaystyle \frac{1}{z^2}}}.`$ (38)
$`r_V^{(1)}`$ can be found in Eq. (18).

The non-Abelian contribution can be written in the form

$`r_V^\mathrm{A}(s,M_Q)`$ $`=`$ $`r_V^{\mathrm{A},pa}+r_V^{\mathrm{A},g}(s,M_Q)|_{\xi =4},`$ (39)
where $`\xi `$ is the QCD gauge parameter, $`r_V^{\mathrm{A},g}`$ is the result from the gluonic double-bubble diagram given below and $`r_V^{\mathrm{A},pa}`$ is again given in terms of a Padé approximant:

$`r_V^{\mathrm{A},pa}`$ $`=`$ $`C_\mathrm{A}C_\mathrm{F}{\displaystyle \frac{3}{4\pi }}\text{Im}\{i\pi {\displaystyle \frac{3327z6z^2}{8z^2\sqrt{11/z}}}+(1+\stackrel{~}{\omega })^2`$ (40)
$`\times [53.56935816521245+23.42800818317322\stackrel{~}{\omega }28.58891755403711\stackrel{~}{\omega }^2`$
$`10.11103466811403\stackrel{~}{\omega }^3+0.4266396586616486\stackrel{~}{\omega }^4]`$
$`\times [19.192245895101771.767625499084376\stackrel{~}{\omega }13.2990319119615\stackrel{~}{\omega }^2`$
$`+1.00664474828887\stackrel{~}{\omega }^3+\stackrel{~}{\omega }^4]^1\}.`$
Following Ref. , the exact results of the imaginary part of the fermionic and gluonic double-bubble diagram (see Fig. 3$`(d)`$ and $`(e)`$), $`r_V^\mathrm{L}`$ and $`r_V^{\mathrm{A},g}`$, respectively, may be written as ($`x\{\mathrm{L};\mathrm{A},g\}`$)

$`r_V^x(s,M_Q)`$ $`=`$ $`C_x\left[{\displaystyle \frac{1}{3}}\left(R_{\mathrm{}}^x\mathrm{ln}{\displaystyle \frac{\mu ^2}{s}}R_0^x\right){\displaystyle \frac{r_V^{(1)}(s,M_Q)}{C_\mathrm{F}}}+R_{\mathrm{}}^x\delta ^{(2)}\right],`$ (41)
with

$`C_\mathrm{L}`$ $`=`$ $`C_\mathrm{F}T,`$
$`C_{\mathrm{A},g}`$ $`=`$ $`C_\mathrm{A}C_\mathrm{F},`$
$`R_{\mathrm{}}^\mathrm{L}`$ $`=`$ $`1,`$
$`R_0^\mathrm{L}`$ $`=`$ $`{\displaystyle \frac{5}{3}}+\mathrm{ln}4,`$
$`R_{\mathrm{}}^{\mathrm{A},g}`$ $`=`$ $`{\displaystyle \frac{5}{4}}{\displaystyle \frac{3}{8}}\xi ,`$
$`R_0^{\mathrm{A},g}`$ $`=`$ $`{\displaystyle \frac{31}{12}}{\displaystyle \frac{3}{4}}\xi +{\displaystyle \frac{3}{16}}\xi ^2+\left({\displaystyle \frac{5}{4}}{\displaystyle \frac{3}{8}}\xi \right)\mathrm{ln}4.`$ (42)
$`r_V^{(1)}`$ is defined in (18), whereas $`\delta ^{(2)}`$ was originally calculated in Ref.  and reads:

$`\delta ^{(2)}={\displaystyle \frac{\left(3v_{Q}^{}{}_{}{}^{2}\right)\left(1+v_{Q}^{}{}_{}{}^{2}\right)}{6}}`$ (43)
$`\times \{\text{Li}_3(p)2\text{Li}_3(1p)3\text{Li}_3(p^2)4\text{Li}_3\left({\displaystyle \frac{p}{1+p}}\right)5\text{Li}_3(1p^2)+{\displaystyle \frac{11}{2}}\zeta _3`$
$`+\text{Li}_2(p)\mathrm{ln}\left({\displaystyle \frac{4\left(1v_{Q}^{}{}_{}{}^{2}\right)}{v_{Q}^{}{}_{}{}^{4}}}\right)+2\text{Li}_2(p^2)\mathrm{ln}\left({\displaystyle \frac{1v_{Q}^{}{}_{}{}^{2}}{2v_{Q}^{}{}_{}{}^{2}}}\right)+2\zeta _2\left[\mathrm{ln}p\mathrm{ln}\left({\displaystyle \frac{1v_{Q}^{}{}_{}{}^{2}}{4v_Q}}\right)\right]`$
$`{\displaystyle \frac{1}{6}}\mathrm{ln}\left({\displaystyle \frac{1+v_Q}{2}}\right)\left[36\mathrm{ln}2\mathrm{ln}p44\mathrm{ln}^2p+49\mathrm{ln}p\mathrm{ln}\left({\displaystyle \frac{1v_{Q}^{}{}_{}{}^{2}}{4}}\right)+\mathrm{ln}^2\left({\displaystyle \frac{1v_{Q}^{}{}_{}{}^{2}}{4}}\right)\right]`$
$`{\displaystyle \frac{1}{2}}\mathrm{ln}p\mathrm{ln}v_Q[36\mathrm{ln}2+21\mathrm{ln}p+16\mathrm{ln}v_Q22\mathrm{ln}(1v_{Q}^{}{}_{}{}^{2})]\}`$
$`+{\displaystyle \frac{1}{24}}\{(156v_{Q}^{}{}_{}{}^{2}v_{Q}^{}{}_{}{}^{4})(\text{Li}_2(p)+\text{Li}_2(p^2))+3(722v_{Q}^{}{}_{}{}^{2}+7v_{Q}^{}{}_{}{}^{4})\text{Li}_2(p)`$
$`(1v_Q)(5145v_Q27v_{Q}^{}{}_{}{}^{2}+5v_{Q}^{}{}_{}{}^{3})\zeta _2`$
$`+{\displaystyle \frac{\left(1+v_Q\right)\left(9+33v_Q9v_{Q}^{}{}_{}{}^{2}15v_{Q}^{}{}_{}{}^{3}+4v_{Q}^{}{}_{}{}^{4}\right)}{v_Q}}\mathrm{ln}^2p`$
$`+[(33+22v_{Q}^{}{}_{}{}^{2}7v_{Q}^{}{}_{}{}^{4})\mathrm{ln}210(3v_{Q}^{}{}_{}{}^{2})(1+v_{Q}^{}{}_{}{}^{2})\mathrm{ln}v_Q`$
$`(1522v_{Q}^{}{}_{}{}^{2}+3v_{Q}^{}{}_{}{}^{4})\mathrm{ln}\left({\displaystyle \frac{1v_{Q}^{}{}_{}{}^{2}}{4v_{Q}^{}{}_{}{}^{2}}}\right)]\mathrm{ln}p`$
$`+2v_Q(3v_{Q}^{}{}_{}{}^{2})\mathrm{ln}\left({\displaystyle \frac{4\left(1v_{Q}^{}{}_{}{}^{2}\right)}{v_{Q}^{}{}_{}{}^{4}}}\right)\left[\mathrm{ln}v_Q3\mathrm{ln}\left({\displaystyle \frac{1v_{Q}^{}{}_{}{}^{2}}{4v_Q}}\right)\right]`$
$`+{\displaystyle \frac{23796v_Q+62v_{Q}^{}{}_{}{}^{2}+32v_{Q}^{}{}_{}{}^{3}59v_{Q}^{}{}_{}{}^{4}}{4}}\mathrm{ln}p16v_Q(3v_{Q}^{}{}_{}{}^{2})\mathrm{ln}\left({\displaystyle \frac{1+v_Q}{4}}\right)`$
$`2v_Q(3917v_{Q}^{}{}_{}{}^{2})\mathrm{ln}\left({\displaystyle \frac{1v_{Q}^{}{}_{}{}^{2}}{2v_{Q}^{}{}_{}{}^{2}}}\right){\displaystyle \frac{v_Q\left(7529v_{Q}^{}{}_{}{}^{2}\right)}{2}}\},`$
with $`\zeta _2=\pi ^2/61.64493`$, $`\zeta _31.20206`$, and

$`p`$ $`=`$ $`{\displaystyle \frac{1v_Q}{1+v_Q}}.`$
The double-bubble contribution $`r_V^\mathrm{T}`$ can be written in the form

$`r_V^\mathrm{T}(s,M_Q)`$ $`=`$ $`TC_\mathrm{F}\left(\rho ^V(M_Q^2,M_Q^2,s)+\rho ^R(M_Q^2,M_Q^2,s)+{\displaystyle \frac{1}{3}}\mathrm{ln}{\displaystyle \frac{M_Q^2}{\mu ^2}}{\displaystyle \frac{r_V^{(1)}}{3C_\mathrm{F}}}\right),`$ (44)
where

$`\rho ^V(M_Q^2,M_Q^2,s)`$ $`=`$ $`{\displaystyle \frac{1}{6}}[{\displaystyle \frac{3+10v_Q^25v_Q^4}{24}}\mathrm{ln}^3p+{\displaystyle \frac{3+40v_Q^2+16v_Q^415v_Q^6}{12v_Q^3}}\mathrm{ln}^2p`$ (45)
$`+\left({\displaystyle \frac{18+234v_Q^2+167v_Q^4118v_Q^6}{18v_Q^2}}+{\displaystyle \frac{310v_Q^2+5v_Q^4}{2}}\zeta _2\right)\mathrm{ln}p`$
$`+{\displaystyle \frac{9+510v_Q^2118v_Q^4}{9v_Q}}+v_Q(27+5v_Q^2)\zeta _2],`$
and

$$\rho ^R(M_Q^2,M_Q^2,s)=\frac{1}{3}_{4M_Q^2/s}^{(12M_Q/\sqrt{s})^2}dy_{4M_Q^2/s}^{(1\sqrt{y})^2}\frac{\mathrm{d}z}{z}\left(1+\frac{2M_Q^2}{sz}\right)\sqrt{1\frac{4M_Q^2}{sz}}(y,z),$$
(46)
with

$`(y,z)`$ $`=`$ $`{\displaystyle \frac{\frac{8M_Q^4}{s^2}+\frac{4M_Q^2}{s}(1y+z)(1y+z)^22(1+z)y}{1y+z}}`$
$`\times \mathrm{ln}{\displaystyle \frac{1y+z\sqrt{1\frac{4M_Q^2}{sy}}\mathrm{\Lambda }^{1/2}(1,y,z)}{1y+z+\sqrt{1\frac{4M_Q^2}{sy}}\mathrm{\Lambda }^{1/2}(1,y,z)}}`$
$`\sqrt{1{\displaystyle \frac{4M_Q^2}{sy}}}\mathrm{\Lambda }^{1/2}(1,y,z)\left[1+{\displaystyle \frac{\frac{16M_Q^4}{s^2}+\frac{8M_Q^2}{s}+4\left(1+\frac{2M_Q^2}{s}\right)z}{(1y+z)^2\left(1\frac{4M_Q^2}{sy}\right)\mathrm{\Lambda }(1,y,z)}}\right],`$
$`\mathrm{\Lambda }(1,y,z)`$ $`=`$ $`1+y^2+z^22(y+z+yz).`$ (47)
Note that $`\rho ^R(M_Q^2,M_Q^2,s)`$ vanishes for $`s(4M_Q)^2`$.

The numerical integration of Eq. (46) would be a bottle neck of rhad. But it turns out that, for $`s>(4M_Q)^2`$, $`r_V^\mathrm{T}(s,M_Q)`$ is extremely well approximated by the high energy expansion of Ref. , and it drastically shortens the running time of rhad. Thus we use this expression as the default in rhad for $`s>(4M_Q)^2`$, whereas for $`s(4M_Q)^2`$ we can savely use Eq. (44) since the double integral does not contribute in this case. For completeness, we list the expansion for $`r_V^\mathrm{T}`$ up to $`M_Q^{12}/s^6`$ :

$$\begin{array}{cc}\hfill r_V^\mathrm{T}(s,M_Q)& =C_\mathrm{F}T\{\frac{11}{8}+\zeta _3+\frac{1}{4}L_{s\mu }+\frac{M_Q^2}{s}[\frac{13}{2}+3L_{s\mu }]\hfill \\ & +\left(\frac{M_Q^2}{s}\right)^2\left[\frac{2}{3}+18\zeta _210\zeta _3+12L_{ms}3L_{ms}^2+\left(\frac{5}{2}6L_{ms}\right)L_{s\mu }\right]\hfill \\ & +\left(\frac{M_Q^2}{s}\right)^3\left[\frac{4}{3}+\frac{352}{9}\zeta _2+\frac{350}{9}L_{ms}\frac{76}{9}L_{ms}^2+\left(\frac{188}{27}\frac{116}{9}L_{ms}\right)L_{s\mu }\right]\hfill \\ & +\left(\frac{M_Q^2}{s}\right)^4[\frac{20233}{288}+\frac{533}{6}\zeta _2+\frac{1673}{72}L_{ms}\frac{25}{4}L_{ms}^2\hfill \\ & \text{ }+(\frac{983}{36}\frac{203}{6}L_{ms})L_{s\mu }]\hfill \\ & +\left(\frac{M_Q^2}{s}\right)^5[\frac{54559}{6750}+\frac{3592}{45}\zeta _2\frac{1157}{75}L_{ms}+\frac{4328}{45}L_{ms}^2\hfill \\ & \text{ }+(\frac{61699}{675}\frac{4676}{45}L_{ms})L_{s\mu }]\hfill \\ & +\left(\frac{M_Q^2}{s}\right)^6[\frac{9214697}{6480}1105\zeta _2+\frac{346981}{540}L_{ms}+\frac{18937}{18}L_{ms}^2\hfill \\ & \text{ }+(\frac{84743}{270}\frac{3064}{9}L_{ms})L_{s\mu }]\}+\mathrm{},\hfill \end{array}$$
(48)
where $`L_{ms}=\mathrm{ln}(M_Q^2/s)`$, $`L_{s\mu }=\mathrm{ln}(s/\mu ^2)`$. The ellipse denotes uncalculated higher order terms in $`M_Q^2/s`$.

The double-bubble contribution with zero outer mass and inner mass $`M_2`$ is given by

$`r_V^{\mathrm{db}}(s,0,M_2)=C_\mathrm{F}T\left[\varrho ^V(M_2^2,s)+\theta (s4M_2^2)\left(\varrho ^R(M_2^2,s)+{\displaystyle \frac{1}{4}}\mathrm{ln}{\displaystyle \frac{M_2^2}{\mu ^2}}\right)\right],`$ (49)
with

$`\varrho ^R(M_2^2,s)`$ $`=`$ $`{\displaystyle \frac{4}{3}}(16x^2)[{\displaystyle \frac{1}{2}}\mathrm{Li}_3({\displaystyle \frac{1w}{2}}){\displaystyle \frac{1}{2}}\mathrm{Li}_3({\displaystyle \frac{1+w}{2}})`$ (50)
$`+\mathrm{Li}_3({\displaystyle \frac{1+w}{1+a}})\mathrm{Li}_3({\displaystyle \frac{1w}{1a}})+\mathrm{Li}_3({\displaystyle \frac{1+w}{1a}})\mathrm{Li}_3({\displaystyle \frac{1w}{1+a}})`$
$`+{\displaystyle \frac{1}{2}}\mathrm{ln}({\displaystyle \frac{1+w}{1w}})\{\zeta _2{\displaystyle \frac{1}{12}}\mathrm{ln}^2({\displaystyle \frac{1+w}{1w}})+{\displaystyle \frac{1}{2}}\mathrm{ln}^2({\displaystyle \frac{a1}{a+1}})`$
$`{\displaystyle \frac{1}{2}}\mathrm{ln}({\displaystyle \frac{1+w}{2}})\mathrm{ln}({\displaystyle \frac{1w}{2}})\}]`$
$`+{\displaystyle \frac{1}{9}}a(19+46x)(\mathrm{Li}_2({\displaystyle \frac{1+w}{1+a}})+\mathrm{Li}_2({\displaystyle \frac{1w}{1a}})\mathrm{Li}_2({\displaystyle \frac{1+w}{1a}})\mathrm{Li}_2({\displaystyle \frac{1w}{1+a}})`$
$`+\mathrm{ln}({\displaystyle \frac{a1}{a+1}})\mathrm{ln}({\displaystyle \frac{1+w}{1w}}))`$
$`+\mathrm{\hspace{0.17em}4}\left({\displaystyle \frac{19}{72}}+x+x^2\right)\left(\mathrm{Li}_2({\displaystyle \frac{1+w}{1w}})\mathrm{Li}_2({\displaystyle \frac{1w}{1+w}})\mathrm{ln}x\mathrm{ln}({\displaystyle \frac{1+w}{1w}})\right)`$
$`+\mathrm{\hspace{0.17em}7}\left({\displaystyle \frac{73}{189}}+{\displaystyle \frac{74}{63}}x+x^2\right)\mathrm{ln}({\displaystyle \frac{1+w}{1w}}){\displaystyle \frac{1}{3}}\left({\displaystyle \frac{2123}{108}}+{\displaystyle \frac{2489}{54}}x\right)w,`$
$`\varrho ^V(M_2^2,s)`$ $`=`$ $`{\displaystyle \frac{2}{3}}\left(16x^2\right)\left(\mathrm{Li}_3(A^2)\zeta _32\zeta _2\mathrm{ln}A+{\displaystyle \frac{2}{3}}\mathrm{ln}^3A\right)`$ (51)
$`+{\displaystyle \frac{1}{9}}\left(19+46x\right)\sqrt{1+4x}\left(\mathrm{Li}_2(A^2)\zeta _2+\mathrm{ln}^2A\right)`$
$`+{\displaystyle \frac{5}{36}}\left({\displaystyle \frac{53}{3}}+44x\right)\mathrm{ln}x+{\displaystyle \frac{3355}{648}}+{\displaystyle \frac{119}{9}}x,`$
where

$`x={\displaystyle \frac{M_2^2}{s}},a=\sqrt{1+4x},w=\sqrt{14x},A={\displaystyle \frac{\sqrt{1+4x}1}{\sqrt{4x}}}.`$
In the limit $`sM_2^2`$, the result for $`r_V^{\mathrm{db}}(s,0,M_2)`$ reads

$`r_V^{\mathrm{db}}(s,0,M_2)`$ $`=`$ $`C_\mathrm{F}T\left({\displaystyle \frac{11}{8}}+\zeta _3+{\displaystyle \frac{1}{4}}\mathrm{ln}{\displaystyle \frac{s}{\mu ^2}}+𝒪\left({\displaystyle \frac{M_2^4}{s^2}}\right)\right).`$ (52)
Furthermore we need $`r_V^{\mathrm{db}}(s,M_1,0)`$ up to the quadratic mass terms which enters the evaluation of $`R_c^{(2)}`$ in the limit $`s>s_b^{\mathrm{thr}}`$. It is given by

$`r_V^{\mathrm{db}}(s,M_1,0)`$ $`=`$ $`C_\mathrm{F}T\left({\displaystyle \frac{11}{8}}+\zeta _3+{\displaystyle \frac{1}{4}}\mathrm{ln}{\displaystyle \frac{s}{\mu ^2}}+{\displaystyle \frac{M_1^2}{s}}\left({\displaystyle \frac{13}{2}}+3\mathrm{ln}{\displaystyle \frac{s}{\mu ^2}}\right)+𝒪\left({\displaystyle \frac{M_1^4}{s^2}}\right)\right).`$
In order to obtain the result for the double-bubble contribution in the limit $`M_2^2s`$ one has to perform the decoupling of $`M_2`$ from the coupling constant using Eq. (9) and the one-loop approximation of Eq. (65). Taking into account the full $`M_1`$ dependence the result can be cast in the form

$`r_V^{\mathrm{db},\mathrm{h}}(s,M_1,M_2)`$ $`=`$ $`C_\mathrm{F}T{\displaystyle \frac{s}{45M_2^2}}[(1+{\displaystyle \frac{12M_1^4}{s^2}}+{\displaystyle \frac{16M_1^6}{s^3}})\mathrm{ln}{\displaystyle \frac{1+\sqrt{1\frac{4M_1^2}{s}}}{1\sqrt{1\frac{4M_1^2}{s}}}}`$ (54)
$`+\left({\displaystyle \frac{22}{5}}+{\displaystyle \frac{79M_1^2}{5s}}+{\displaystyle \frac{8M_1^4}{s^2}}\right)\sqrt{1{\displaystyle \frac{4M_1^2}{s}}}`$
$`+(1+{\displaystyle \frac{2M_1^2}{s}})\mathrm{ln}{\displaystyle \frac{M_2^2}{M_1^2}}\sqrt{1{\displaystyle \frac{4M_1^2}{s}}}].`$
The expansion of $`C_\mathrm{F}T\varrho ^V`$ in the limit $`M_2^2s`$ agrees with $`r_V^{\mathrm{db},\mathrm{h}}`$ for $`M_1=0`$ and reads

$`r_V^{\mathrm{db},\mathrm{h}}(s,0,M_2)`$ $`=`$ $`C_\mathrm{F}T{\displaystyle \frac{3s}{8M_2^2}}\left({\displaystyle \frac{176}{675}}+{\displaystyle \frac{8}{135}}\mathrm{ln}{\displaystyle \frac{M_2^2}{s}}\right).`$ (55)
The massless limits of the above expressions have been given in Eq. (22).

### E.2 Analytical results for $`R_Q^{(3)}(s)`$

This section lists analytical expressions for the three-loop functions defined in Sect. 6. As already mentioned, only the high-energy expansion is known, including quartic mass terms. We write these functions as follows:

$$\begin{array}{cc}\hfill r_V^{(3)}& =r_0^{(3)}+\frac{M_Q^2}{s}r_{V,2}^{(3)}+\left(\frac{M_Q^2}{s}\right)^2r_{V,4}^{(3)},\delta r_0^{(3)}=\frac{M_Q^2}{s}r_{0,2}^{(3)}+\left(\frac{M_Q^2}{s}\right)^2r_{0,4}^{(3)}.\hfill \end{array}$$
(56)
The results for the massless terms have been obtained in Ref. , the quadratic mass terms in Ref. , and the quartic terms in Ref. . Adopting the on-shell scheme for the mass $`M_Q`$, we get, for the contributions where the massive quark couples to the external current:

$$\begin{array}{cc}\hfill r_0^{(3)}& =\frac{87029}{288}\frac{121}{8}\zeta _2\frac{1103}{4}\zeta _3+\frac{275}{6}\zeta _5+L_{s\mu }\left(\frac{4321}{48}+\frac{121}{2}\zeta _3\right)+\frac{121}{16}L_{s\mu }^2\hfill \\ & +n_f\left[\frac{7847}{216}+\frac{11}{6}\zeta _2+\frac{262}{9}\zeta _3\frac{25}{9}\zeta _5+L_{s\mu }\left(\frac{785}{72}\frac{22}{3}\zeta _3\right)\frac{11}{12}L_{s\mu }^2\right]\hfill \\ & +n_f^2\left[\frac{151}{162}\frac{1}{18}\zeta _2\frac{19}{27}\zeta _3+L_{s\mu }\left(\frac{11}{36}+\frac{2}{9}\zeta _3\right)+\frac{1}{36}L_{s\mu }^2\right],\hfill \end{array}$$
(57)
$$\begin{array}{cc}\hfill r_{V,2}^{(3)}& =\frac{22351}{12}\frac{967}{2}\zeta _216\zeta _2\mathrm{ln}2+\frac{502}{3}\zeta _3\frac{5225}{6}\zeta _5+378L_{ms}9L_{ms}^2\hfill \\ & +L_{s\mu }(\frac{2385}{4}132L_{ms})+\frac{363}{4}L_{s\mu }^2+n_f[\frac{8429}{54}+42\zeta _2\frac{466}{27}\zeta _3+\frac{1045}{27}\zeta _5\hfill \\ & \frac{52}{3}L_{ms}+2L_{ms}^2+L_{s\mu }(\frac{389}{6}+8L_{ms})11L_{s\mu }^2]\hfill \\ & +n_f^2\left(\frac{125}{54}\frac{2}{3}\zeta _2\frac{13}{9}L_{s\mu }+\frac{1}{3}L_{s\mu }^2\right),\hfill \end{array}$$
(58)
$$\begin{array}{cc}\hfill r_{V,4}^{(3)}& =\frac{91015}{108}\frac{76}{9}\mathrm{ln}^42+\frac{2564287}{540}\zeta _2\frac{4568}{9}\zeta _2\mathrm{ln}2\hfill \\ & \frac{128}{3}\zeta _2\mathrm{ln}^22+\frac{56257}{18}\zeta _3\frac{1439}{3}\zeta _2\zeta _3\frac{1565}{6}\zeta _4\frac{3770}{3}\zeta _5\frac{608}{3}a_4\hfill \\ & +L_{ms}\left(\frac{5536}{3}+564\zeta _224\zeta _2\mathrm{ln}2+416\zeta _3\right)\frac{591}{4}L_{ms}^2+\frac{15}{2}L_{ms}^3\hfill \\ & +L_{s\mu }\left(\frac{5297}{12}1199\zeta _2572\zeta _388\zeta _2\mathrm{ln}2+\frac{4033}{4}L_{ms}+\frac{165}{2}L_{ms}^2\right)\hfill \\ & +L_{s\mu }^2\left(\frac{605}{8}\frac{363}{2}L_{ms}\right)\hfill \\ & +n_f[\frac{21011}{216}+\frac{8}{27}\mathrm{ln}^42\frac{3544}{9}\zeta _2\frac{176}{9}\zeta _2\mathrm{ln}2+\frac{32}{9}\zeta _2\mathrm{ln}^22\frac{2323}{9}\zeta _3\hfill \\ & +\frac{700}{9}\zeta _4+\frac{440}{9}\zeta _5+\frac{64}{9}a_4+L_{ms}\left(\frac{2419}{12}+\frac{44}{3}\zeta _2+\frac{16}{3}\zeta _2\mathrm{ln}2+\frac{28}{3}\zeta _3\right)\hfill \\ & \frac{157}{6}L_{ms}^2\frac{2}{3}L_{ms}^3+L_{s\mu }(\frac{1879}{36}+\frac{416}{3}\zeta _2+\frac{148}{3}\zeta _3+\frac{16}{3}\zeta _2\mathrm{ln}2\hfill \\ & \frac{361}{3}L_{ms}+6L_{ms}^2)+L_{s\mu }^2(\frac{55}{6}+22L_{ms})]\hfill \\ & +n_f^2[\frac{35}{18}+\frac{25}{3}\zeta _2+\frac{112}{27}\zeta _3+L_{ms}(\frac{94}{27}\frac{8}{3}\zeta _2)+\frac{13}{9}L_{ms}^2\frac{2}{9}L_{ms}^3\hfill \\ & +L_{s\mu }(\frac{35}{27}4\zeta _2\frac{8}{9}\zeta _3+3L_{ms}\frac{2}{3}L_{ms}^2)+L_{s\mu }^2(\frac{5}{18}\frac{2}{3}L_{ms})],\hfill \end{array}$$
(59)
where $`n_f`$ is the number of active flavors, $`a_4=\text{Li}_4(1/2)0.517479`$, $`\zeta _2=\pi ^2/61.64492`$, $`\zeta _31.20206`$, $`\zeta _4=\pi ^4/901.08232`$, $`\zeta _51.03693`$, $`L_{ms}=\mathrm{ln}(M_Q^2/s)`$ and $`L_{s\mu }=\mathrm{ln}(s/\mu ^2)`$. The mass corrections in the case where the massive quark only appears as insertion in gluon lines reads:

$$\begin{array}{cc}\hfill r_{0,2}^{(3)}& =80+60\zeta _3+n_f\left(\frac{32}{9}\frac{8}{3}\zeta _3\right),\hfill \end{array}$$
(60)
$$\begin{array}{cc}\hfill r_{0,4}^{(3)}& =\frac{4217}{48}+15\zeta _2+\frac{139}{3}\zeta _3+\frac{50}{3}\zeta _5+L_{ms}\left(\frac{97}{4}38\zeta _3\right)2L_{ms}^2\hfill \\ & +L_{s\mu }\left(\frac{143}{6}+22\zeta _3+\frac{11}{2}L_{ms}\right)\hfill \\ & +n_f\left[\frac{457}{108}\frac{2}{3}\zeta _2\frac{22}{9}\zeta _3+L_{ms}\left(\frac{13}{18}+\frac{4}{3}\zeta _3\right)+L_{s\mu }\left(\frac{13}{9}\frac{1}{3}L_{ms}\frac{4}{3}\zeta _3\right)\right].\hfill \end{array}$$
(61)
### E.3 Analytical results for $`R_Q^{(4)}(s)`$

With the help of the renormalization group equation the logarithmic contributions of the massless order $`\alpha _s^4`$ result can be reconstructed exactly. Thus, for $`\mu ^2s`$ Eqs. (33) and (34) take the form

$$\begin{array}{cc}\hfill \overline{r}_0^{(4),0\mathrm{L}}(s)& =186+\left(\frac{520175}{192}+\frac{3993}{32}\zeta _2+\frac{38643}{16}\zeta _3\frac{3025}{8}\zeta _5\right)L_{s\mu }\hfill \\ & +\left(\frac{49775}{128}\frac{3993}{16}\zeta _3\right)L_{s\mu }^2\frac{1331}{64}L_{s\mu }^3,\hfill \\ \hfill \overline{r}_0^{(4),1\mathrm{L}}(s)& =21.3+\left(\frac{188521}{384}\frac{363}{16}\zeta _2\frac{9695}{24}\zeta _3+\frac{275}{6}\zeta _5\right)L_{s\mu }\hfill \\ & +\left(\frac{2263}{32}+\frac{363}{8}\zeta _3\right)L_{s\mu }^2+\frac{121}{32}L_{s\mu }^3,\hfill \\ \hfill r_0^{(4),2\mathrm{L}}(s)& =\frac{1045381}{15552}\frac{593}{432}\pi ^2\frac{40655}{864}\zeta _3+\frac{11}{12}\pi ^2\zeta _3+\frac{5}{6}\zeta _3^2\frac{260}{27}\zeta _5\hfill \\ & +\left(\frac{94693}{3456}+\frac{11}{8}\zeta _2+\frac{257}{12}\zeta _3\frac{25}{18}\zeta _5\right)L_{s\mu }+\left(\frac{593}{144}\frac{11}{4}\zeta _3\right)L_{s\mu }^2\frac{11}{48}L_{s\mu }^3,\hfill \\ \hfill r_0^{(4),3\mathrm{L}}(s)& =\frac{6131}{5832}+\frac{11}{432}\pi ^2+\frac{203}{324}\zeta _3\frac{1}{54}\pi ^2\zeta _3+\frac{5}{18}\zeta _5\hfill \\ & +\left(\frac{151}{324}\frac{1}{36}\zeta _2\frac{19}{54}\zeta _3\right)L_{s\mu }+\left(\frac{11}{144}+\frac{1}{18}\zeta _3\right)L_{s\mu }^2+\frac{1}{216}L_{s\mu }^3,\hfill \\ \hfill \overline{r}_{V,2}^{(4),0\mathrm{L}}(s)& =7110\frac{9855433}{864}\frac{535211}{270}\zeta _2+\frac{3194}{9}\zeta _2\mathrm{ln}2+\frac{128}{3}\zeta _2\mathrm{ln}^22+\frac{76}{9}\mathrm{ln}^42\frac{9655}{18}\zeta _3\hfill \\ & +\frac{1439}{3}\zeta _2\zeta _3+\frac{1565}{6}\zeta _4+\frac{18925}{9}\zeta _5+\frac{608}{3}a_4\hfill \\ & +\left(\frac{178849}{24}659\zeta _2+56\zeta _2\mathrm{ln}2+\frac{938}{3}\zeta _3\frac{5225}{3}\zeta _5\right)L_{ms}\frac{2313}{8}L_{ms}^2+\frac{21}{2}L_{ms}^3\hfill \\ & +(\frac{538849}{32}\frac{6849}{2}L_{ms}+\frac{297}{4}L_{ms}^2+\frac{31911}{8}\zeta _2+132\zeta _2\mathrm{ln}2\frac{2761}{2}\zeta _3\hfill \\ & +\frac{57475}{8}\zeta _5)L_{s\mu }+(\frac{85437}{32}+\frac{1089}{2}L_{ms})L_{s\mu }^2\frac{3993}{16}L_{s\mu }^3,\hfill \\ \hfill \overline{r}_{V,2}^{(4),1\mathrm{L}}(s)& =1430+\frac{604967}{648}+\frac{689}{3}\zeta _2+\frac{76}{3}\zeta _2\mathrm{ln}2\frac{32}{9}\zeta _2\mathrm{ln}^22\frac{8}{27}\mathrm{ln}^42+\frac{10622}{81}\zeta _3\hfill \\ & \frac{610}{9}\zeta _4\frac{8360}{81}\zeta _5\frac{64}{9}a_4+(\frac{35899}{54}+\frac{64}{3}\zeta _2\frac{16}{3}\zeta _2\mathrm{ln}2\frac{1436}{27}\zeta _3\hfill \\ & +\frac{2090}{27}\zeta _5)L_{ms}+\frac{133}{2}L_{ms}^2\frac{10}{3}L_{ms}^3+(\frac{713191}{288}\frac{2353}{4}\zeta _28\zeta _2\mathrm{ln}2+\frac{4069}{18}\zeta _3\hfill \\ & \frac{13585}{18}\zeta _5+370L_{ms}21L_{ms}^2)L_{s\mu }+(\frac{3643}{8}66L_{ms})L_{s\mu }^2+\frac{363}{8}L_{s\mu }^3,\hfill \\ \hfill r_{V,2}^{(4),2\mathrm{L}}(s)& =\frac{114253}{648}\frac{415}{6}\zeta _2\frac{995}{324}\zeta _3+\frac{53}{9}\zeta _3^2\frac{3610}{81}\zeta _5+\left(\frac{199}{18}+\frac{4}{3}\zeta _2\right)L_{ms}\frac{13}{6}L_{ms}^2\hfill \\ & +\frac{2}{9}L_{ms}^3+\left(\frac{90835}{864}+\frac{53}{2}\zeta _2\frac{233}{27}\zeta _3+\frac{1045}{54}\zeta _5\frac{26}{3}L_{ms}+L_{ms}^2\right)L_{s\mu }\hfill \\ & +\left(\frac{95}{4}+2L_{ms}\right)L_{s\mu }^2\frac{11}{4}L_{s\mu }^3,\hfill \\ \hfill r_{V,2}^{(4),3\mathrm{L}}(s)& =\frac{2705}{1944}+\frac{13}{18}\zeta _2+\left(\frac{125}{108}\frac{1}{3}\zeta _2\right)L_{s\mu }\frac{13}{36}L_{s\mu }^2+\frac{1}{18}L_{s\mu }^3,\hfill \end{array}$$
(62)
where the first number in $`\overline{r}_0^{(4),0\mathrm{L}}`$, $`\overline{r}_0^{(4),1\mathrm{L}}`$, $`\overline{r}_{V,2}^{(4),0\mathrm{L}}`$, and $`\overline{r}_{V,2}^{(4),1\mathrm{L}}`$ corresponds to the estimates obtained with the help of the results of Refs. .

### E.4 Renormalization group coefficients

This section collects the formulas used for the evaluation of the running coupling constant and the $`\overline{\text{MS}}`$ quark masses. It also contains the conversion formulas for transforming $`\overline{\text{MS}}`$ quark masses to their on-shell values.

The coefficients of the $`\beta `$-function defined in Eq. (8) read :

$`\beta _0^{(n_f)}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left[11{\displaystyle \frac{2}{3}}n_f\right],`$
$`\beta _1^{(n_f)}`$ $`=`$ $`{\displaystyle \frac{1}{16}}\left[102{\displaystyle \frac{38}{3}}n_f\right],`$
$`\beta _2^{(n_f)}`$ $`=`$ $`{\displaystyle \frac{1}{64}}\left[{\displaystyle \frac{2857}{2}}{\displaystyle \frac{5033}{18}}n_f+{\displaystyle \frac{325}{54}}n_f^2\right],`$
$`\beta _3^{(n_f)}`$ $`=`$ $`{\displaystyle \frac{1}{256}}[({\displaystyle \frac{149753}{6}}+3564\zeta _3)({\displaystyle \frac{1078361}{162}}+{\displaystyle \frac{6508}{27}}\zeta _3)n_f`$ (63)
$`+({\displaystyle \frac{50065}{162}}+{\displaystyle \frac{6472}{81}}\zeta _3)n_f^2+{\displaystyle \frac{1093}{729}}n_f^3].`$
The coefficients of the $`\gamma _m`$-function of Eq. (10) are :

$`\gamma _{m,0}^{(n_f)}`$ $`=`$ $`1,`$
$`\gamma _{m,1}^{(n_f)}`$ $`=`$ $`{\displaystyle \frac{1}{16}}\left[{\displaystyle \frac{202}{3}}{\displaystyle \frac{20}{9}}n_f\right],`$
$`\gamma _{m,2}^{(n_f)}`$ $`=`$ $`{\displaystyle \frac{1}{64}}\left[1249+\left({\displaystyle \frac{2216}{27}}{\displaystyle \frac{160}{3}}\zeta _3\right)n_f{\displaystyle \frac{140}{81}}n_f^2\right],`$
$`\gamma _{m,3}^{(n_f)}`$ $`=`$ $`{\displaystyle \frac{1}{256}}[{\displaystyle \frac{4603055}{162}}+{\displaystyle \frac{135680}{27}}\zeta _38800\zeta _5`$ (64)
$`+\left({\displaystyle \frac{91723}{27}}{\displaystyle \frac{34192}{9}}\zeta _3+880\zeta _4+{\displaystyle \frac{18400}{9}}\zeta _5\right)n_f`$
$`+({\displaystyle \frac{5242}{243}}+{\displaystyle \frac{800}{9}}\zeta _3{\displaystyle \frac{160}{3}}\zeta _4)n_f^2+({\displaystyle \frac{332}{243}}+{\displaystyle \frac{64}{27}}\zeta _3)n_f^3].`$
The decoupling coefficient for $`\alpha _s`$ of Eq. (9) is :

$`\left(\zeta _g\right)^2`$ $`=`$ $`1+{\displaystyle \frac{\alpha _s^{(n_f)}(\mu )}{\pi }}\left({\displaystyle \frac{1}{6}}\mathrm{ln}{\displaystyle \frac{\mu ^2}{M_h^2}}\right)+\left({\displaystyle \frac{\alpha _s^{(n_f)}(\mu )}{\pi }}\right)^2\left({\displaystyle \frac{7}{24}}{\displaystyle \frac{19}{24}}\mathrm{ln}{\displaystyle \frac{\mu ^2}{M_h^2}}+{\displaystyle \frac{1}{36}}\mathrm{ln}^2{\displaystyle \frac{\mu ^2}{M_h^2}}\right)`$ (65)
$`+\left({\displaystyle \frac{\alpha _s^{(n_f)}(\mu )}{\pi }}\right)^3[{\displaystyle \frac{58933}{124416}}{\displaystyle \frac{2}{3}}\zeta _2(1+{\displaystyle \frac{1}{3}}\mathrm{ln}2){\displaystyle \frac{80507}{27648}}\zeta _3{\displaystyle \frac{8521}{1728}}\mathrm{ln}{\displaystyle \frac{\mu ^2}{M_h^2}}`$
$`{\displaystyle \frac{131}{576}}\mathrm{ln}^2{\displaystyle \frac{\mu ^2}{M_h^2}}{\displaystyle \frac{1}{216}}\mathrm{ln}^3{\displaystyle \frac{\mu ^2}{M_h^2}}+n_l({\displaystyle \frac{2479}{31104}}+{\displaystyle \frac{\zeta _2}{9}}+{\displaystyle \frac{409}{1728}}\mathrm{ln}{\displaystyle \frac{\mu ^2}{M_h^2}})],`$
and for the quark masses, Eq. (13) :

$`\zeta _m`$ $`=`$ $`1+\left({\displaystyle \frac{\alpha _s^{(n_f)}(\mu )}{\pi }}\right)^2({\displaystyle \frac{89}{432}}{\displaystyle \frac{5}{36}}\mathrm{ln}{\displaystyle \frac{\mu ^2}{M_h^2}}+{\displaystyle \frac{1}{12}}\mathrm{ln}^2{\displaystyle \frac{\mu ^2}{M_h^2}})+\left({\displaystyle \frac{\alpha _s^{(n_f)}(\mu )}{\pi }}\right)^3[{\displaystyle \frac{1871}{2916}}`$ (66)
$`{\displaystyle \frac{407}{864}}\zeta _3+{\displaystyle \frac{5}{4}}\zeta _4{\displaystyle \frac{1}{36}}B_4+\left({\displaystyle \frac{121}{2592}}{\displaystyle \frac{5}{6}}\zeta _3\right)\mathrm{ln}{\displaystyle \frac{\mu ^2}{M_h^2}}+{\displaystyle \frac{319}{432}}\mathrm{ln}^2{\displaystyle \frac{\mu ^2}{M_h^2}}`$
$`+{\displaystyle \frac{29}{216}}\mathrm{ln}^3{\displaystyle \frac{\mu ^2}{M_h^2}}+n_l({\displaystyle \frac{1327}{11664}}{\displaystyle \frac{2}{27}}\zeta _3{\displaystyle \frac{53}{432}}\mathrm{ln}{\displaystyle \frac{\mu ^2}{M_h^2}}{\displaystyle \frac{1}{108}}\mathrm{ln}^3{\displaystyle \frac{\mu ^2}{M_h^2}})],`$
where $`B_41.76280`$, $`n_l=n_f1`$ is the number of light (massless) quarks and $`M_h`$ is the pole mass of the heavy quark. The versions of Eqs. (65) and (66) where the $`\overline{\text{MS}}`$ mass is used as parameter can be found in Refs. .

The conversion from $`\overline{\text{MS}}`$ to on-shell quark masses reads :

$`{\displaystyle \frac{M_Q}{\overline{m}_Q(\mu )}}`$ $`=`$ $`1+{\displaystyle \frac{\alpha _s^{(n_f)}(\mu )}{\pi }}({\displaystyle \frac{4}{3}}+l_{\mu m})+\left({\displaystyle \frac{\alpha _s^{(n_f)}(\mu )}{\pi }}\right)^2[{\displaystyle \frac{307}{32}}+(2+{\displaystyle \frac{2}{3}}\mathrm{ln}2)\zeta _2{\displaystyle \frac{1}{6}}\zeta _3`$ (67)
$`+{\displaystyle \frac{493}{72}}l_{\mu m}+{\displaystyle \frac{43}{24}}l_{\mu m}^2+n_l({\displaystyle \frac{71}{144}}{\displaystyle \frac{1}{3}}\zeta _2{\displaystyle \frac{13}{36}}l_{\mu m}{\displaystyle \frac{1}{12}}l_{\mu m}^2)]`$
$`+\left({\displaystyle \frac{\alpha _s^{(n_f)}(\mu )}{\pi }}\right)^3[{\displaystyle \frac{8481925}{93312}}+\zeta _2({\displaystyle \frac{652841}{6480}}{\displaystyle \frac{1439}{72}}\zeta _3)+{\displaystyle \frac{58}{27}}\zeta _3{\displaystyle \frac{3475}{432}}\zeta _4`$
$`+{\displaystyle \frac{1975}{216}}\zeta _5{\displaystyle \frac{220}{27}}a_4{\displaystyle \frac{575}{27}}\zeta _2\mathrm{ln}2{\displaystyle \frac{44}{27}}\zeta _2\mathrm{ln}^22{\displaystyle \frac{55}{162}}\mathrm{ln}^42`$
$`+l_{\mu m}\left({\displaystyle \frac{177305}{2592}}+{\displaystyle \frac{37}{3}}\zeta _2{\displaystyle \frac{67}{36}}\zeta _3+{\displaystyle \frac{37}{9}}\zeta _2\mathrm{ln}2\right)+{\displaystyle \frac{19315}{864}}l_{\mu m}^2+{\displaystyle \frac{1591}{432}}l_{\mu m}^3`$
$`+n_l({\displaystyle \frac{231847}{23328}}{\displaystyle \frac{991}{108}}\zeta _2{\displaystyle \frac{241}{72}}\zeta _3+{\displaystyle \frac{305}{108}}\zeta _4+{\displaystyle \frac{8}{27}}a_4{\displaystyle \frac{22}{27}}\zeta _2\mathrm{ln}2`$
$`+{\displaystyle \frac{4}{27}}\zeta _2\mathrm{ln}^22+{\displaystyle \frac{1}{81}}\mathrm{ln}^42+l_{\mu m}\left({\displaystyle \frac{10129}{1296}}{\displaystyle \frac{49}{18}}\zeta _2{\displaystyle \frac{7}{9}}\zeta _3{\displaystyle \frac{2}{9}}\zeta _2\mathrm{ln}2\right)`$
$`{\displaystyle \frac{1103}{432}}l_{\mu m}^2{\displaystyle \frac{10}{27}}l_{\mu m}^3)+n_l^2({\displaystyle \frac{2353}{23328}}+{\displaystyle \frac{13}{54}}\zeta _2+{\displaystyle \frac{7}{54}}\zeta _3+l_{\mu m}({\displaystyle \frac{89}{648}}+{\displaystyle \frac{1}{9}}\zeta _2)`$
$`+{\displaystyle \frac{13}{216}}l_{\mu m}^2+{\displaystyle \frac{1}{108}}l_{\mu m}^3)],`$
where $`a_4=\mathrm{Li}_4(1/2)0.517479`$, $`l_{\mu m}=\mathrm{ln}(\mu ^2/\overline{m}_Q^2(\mu ))`$, and $`n_l=n_f1`$. Other variants of Eq. (67) are listed in Ref. .


text:
# On the Superconductivity in the Induced Pairing Model

## 1 The Model

We consider the model of coexisting electron pairs (formed by $`\mathrm{"}d\mathrm{"}`$ electrons) and itinerant $`\mathrm{"}c\mathrm{"}`$ electrons defined by the following effective Hamiltonian

$`H={\displaystyle \underset{𝐤\sigma }{}}(ϵ_𝐤\mu )c_{𝐤\sigma }^{}c_{𝐤\sigma }+2{\displaystyle \underset{i}{}}(\mathrm{\Delta }_0\mu )b_i^{}b_i{\displaystyle \underset{ij}{}}J_{ij}b_i^{}b_j+`$
$`{\displaystyle \underset{𝐤,𝐪}{}}[V_𝐪(𝐤)c_{𝐤+𝐪/\mathrm{𝟐},}^{}c_{𝐤+𝐪/2,}^{}b_𝐪+H.c.],`$ (1)
where $`ϵ_𝐤`$ refers to the band energy of the c-electrons, $`\mathrm{\Delta }_0`$ measures the relative position of the LP level with respect to the bottom of the c-electron band, $`\mu `$ is the chemical potential which ensures that the total number of particles in the system is constant, i.e. $`n=\frac{1}{N}(_{𝐤\sigma }<c_{𝐤\sigma }^{}c_{𝐤\sigma }>+2_i<b_i^{}b_i>)=n_c+2n_B.`$ $`n_c`$ is the concentration of $`c`$-electrons, $`n_B`$ is the number of $`d`$-pairs (hard-core bosons) per site. $`J_{ij}`$ is the pair hopping integral The operators for local pairs $`b_i^{},b_i`$ obey the Pauli spin 1/2 commutation rules. $`V_𝐪(𝐤)`$ describes the coupling between the two subsystems. Furthermore,we set $`J_{ij}=0`$. We will analyze the case $`V_𝐪(𝐤)=V_0(𝐤)=I\varphi _k/\sqrt{N}`$, and neglect its $`q`$ dependence at small $`q`$. The interaction term takes the form of coupling, via the center of mass momenta $`𝐪`$, of the singlet pair of $`c`$-electrons $`B_𝐪^{}`$ and the hard-core boson $`b_𝐪`$:

$`H_1={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{𝐪}{}}I(B_𝐪^{}b_𝐪+b_𝐪^{}B_𝐪).`$ (2)
$`B_𝐪^{}=_𝐤\varphi _kc_{𝐤+𝐪/\mathrm{𝟐},}^{}c_{𝐤+𝐪/\mathrm{𝟐},}^{}`$ denotes the singlet pair creation operator of $`c`$-electrons and $`I`$ is the coupling constant. The pairing symmetry, on a 2D square lattice, is determined by the form of $`\varphi _k`$, which is $`1`$ for on-site pairing ($`s`$), $`\varphi _k=\gamma _k=\mathrm{cos}(k_xa)+\mathrm{cos}(k_ya)`$ for extended $`s`$-wave ($`s^{}`$) and $`\varphi _k=\eta _k=\mathrm{cos}(k_xa)\mathrm{cos}(k_ya)`$ for $`d_{x^2y^2}`$-wave pairing ($`d`$). In general, one can consider a decomposition $`I\varphi _k=g_0+g_s\gamma _k+g_d\eta _k`$, with appropriate coupling parameters for different symmetry channels.

The numerical results presented below are based on the BCS-Mean-Field Approximation (MFA) and the Kosterlitz-Thouless (KT) theory for 2D superfluid. In the BCS-MFA the free energy of the system is calculated to be:

$`F/N={\displaystyle \frac{2}{\beta N}}{\displaystyle \underset{𝐤}{}}\mathrm{ln}\left[2\mathrm{cosh}(\beta E_𝐤/2)\right]{\displaystyle \frac{1}{\beta }}\mathrm{ln}\left[2\mathrm{cosh}(\beta \mathrm{\Delta })\right]+C,`$ (3)
$`C=ϵ_b+\mathrm{\Delta }_0+\mu (n_c+2n_B)2\mu 2I|x_0|\rho _0^x,`$ (4)
where the quasiparticle energy of the $`c`$-electron subsystem is given by
$`E_𝐤=\sqrt{\overline{ϵ}_𝐤^2+\overline{\mathrm{\Delta }}_𝐤^2}`$, and $`\overline{ϵ}_𝐤=ϵ_𝐤\mu `$, $`\overline{\mathrm{\Delta }}_𝐤^2=I^2\varphi _k^2(\rho _0^x)^2`$, $`\mathrm{\Delta }=\sqrt{(\mathrm{\Delta }_0\mu )^2+I^2|x_0|^2}`$, $`\beta =1/k_BT.`$ The $`c`$-electron dispersion is $`ϵ_𝐤=\stackrel{~}{ϵ}_𝐤ϵ_b=2t\left[\mathrm{cos}(k_xa)+\mathrm{cos}(k_ya)\right]4t_2\mathrm{cos}(k_xa)\mathrm{cos}(k_ya)ϵ_b`$, with the nn and nnn hopping parameters $`t`$ and $`t_2`$, respectively, $`ϵ_b=min\stackrel{~}{ϵ}_𝐤`$ . It should be noted that the energy gap in the $`c`$-band is due to nonzero Bose condensate amplitude ($`|<b>|0`$), and well defined Bogolyubov quasiparticles exist in the superconducting phase. The energy spectrum of the system is characterized by $`E_𝐤`$ and $`\mathrm{\Delta }`$. The superconducting order parameters: $`x_0=\frac{1}{N}_k\varphi _k<c_k^{}c_k^{}>`$ and $`\rho _0^x=\frac{1}{2N}_i<b_i^{}+b_i>`$ and the chemical potential $`\mu `$ are given by

$`{\displaystyle \frac{F}{x_0}}=0,{\displaystyle \frac{F}{\rho _0^x}}=0,{\displaystyle \frac{F}{\mu }}=0.`$ (5)
The superfluid density derived within the linear response method and BCS theory is of the form:

$`\rho _s={\displaystyle \frac{1}{2N}}{\displaystyle \underset{k}{}}\left\{\left({\displaystyle \frac{ϵ_𝐤}{k_x}}\right)^2{\displaystyle \frac{f(E_𝐤)}{E_𝐤}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2ϵ_𝐤}{k_x^2}}\left[1{\displaystyle \frac{\overline{ϵ}_𝐤}{E_𝐤}}\mathrm{tanh}\left({\displaystyle \frac{\beta E_𝐤}{2}}\right)\right]\right\}\text{,}`$ (6)
where $`f(E_k)=1/\left[\mathrm{exp}(\beta E_k)+1\right]`$ is the Fermi-Dirac distribution function. In the local limit: $`\lambda ^2=(16\pi e^2/\mathrm{}^2c^2)\rho _s`$, where $`\lambda `$ is the London penetration depth.
Let us also consider the pair propagator for $`c`$-electron subsystem

$`G_2(𝐪,\omega )={\displaystyle \frac{1}{N}}<<B_𝐪|B_𝐪^{}>>_\omega ,`$ (7)
where $`<<B_𝐪|B_𝐪^{}>>_\omega `$ is the time Fourier transform of the Green’s function: $`i\mathrm{\Theta }(tt^{})<[B_𝐪(t),B_𝐪^{}(t^{})]>`$. In the normal state, using the equation of motion technique and Random Phase Approximation one gets:

$`G_2(𝐪,\omega )={\displaystyle \frac{\chi (𝐪,\omega )}{1v_{eff}(\omega )\chi (𝐪,\omega )}},`$ (8)
where

$`\chi (𝐪,\omega )={\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}\varphi _k^2{\displaystyle \frac{1f(\overline{ϵ}_{k+q/2})f(\overline{ϵ}_{kq/2})}{\omega \overline{ϵ}_{k+q/2}\overline{ϵ}_{kq/2}}},`$ (9)
$`v_{eff}(\omega )={\displaystyle \frac{I^2}{\omega 2(\mathrm{\Delta }_0\mu )}}\mathrm{tanh}\left[\beta (\mathrm{\Delta }_0\mu )\right],`$ (10)
$`\chi (𝐪,\omega )`$ is the pair susceptibility and $`v_{eff}(\omega )`$-the effective interaction mediated by LP’s. Introducing the generalized T-matrix $`\mathrm{\Gamma }(𝐪,\omega )`$ via the relation $`G_2=\chi \mathrm{\Gamma }/v_{eff}`$, one has

$`\mathrm{\Gamma }(𝐪,\omega )={\displaystyle \frac{v_{eff}(\omega )}{1v_{eff}(\omega )\chi (𝐪,\omega )}}.`$ (11)
An instability of the normal phase occcurs at $`\mathrm{\Gamma }^1(0,0)=0`$ (the Thouless criterion), and $`T_c`$ is given by

$`1={\displaystyle \frac{I^2}{2(\mathrm{\Delta }_0\mu )}}\mathrm{tanh}\left[\beta _c^{MFA}(\mathrm{\Delta }_0\mu )\right]{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}\varphi _k^2{\displaystyle \frac{\mathrm{tanh}\left(\beta _c^{MFA}\overline{ϵ}_𝐤/2\right)}{2\overline{ϵ}_𝐤}},`$ (12)
in full agreement with the MFA expression which one obtains from Eqs.(5), if $`x_00,\rho _0^x0`$.

The MFA transition temperature ($`T_c^{MFA}`$) at which the gap amplitude vanishes yields an estimation of the $`c`$-electron pair formation temperature . Due to the fluctuation effects the superconducting phase transition will occur at a critical temperature lower than that given by the BCS-MFA theory. In 2D, $`T_c`$ can be derived within the KT theory for superfluid , which describes the transition in terms of vortex-antivortex pair unbinding. We evaluate $`T_c`$ using the KT relation for the universal jump of the superfluid density $`\rho _s`$ at $`T_c`$  :

$$\frac{2}{\pi }k_BT_c=\rho _s(T_c),$$
(13)
where $`\rho _s(T)`$ is given by Eq.(6) and $`x_0(T),\rho _0^x(T),\mu (T)`$ are given by Eqs.(5). Hence, the critical temperature denoted further by $`T_c^{KT}`$ is determined from the set of four self-consistent equations.

## 2 Results and discussion

We have performed a comprehensive analysis of the phase diagrams and superfluid properties of the model (1) for different pairing symmetries . Below we will discuss these results including the supplementary ones.
In Fig.1 we show the ground state diagrams as a function of $`n`$ and $`\mathrm{\Delta }_0/D`$, plotted for $`|I_0|=0`$ ($`I=|I_0|`$) and two fixed values of $`t_2/t`$. A weak intersubsystem coupling $`|I_0|`$ will not change much characteristic lines of the diagrams, while increasing $`|I_0|`$ expands the range of the LP+E regime . Depending on the relative concentration of ”c” and ”d” electrons we distinguish three essentially different physical situations (in the absence of interactions). For $`n2`$ it will be:
(i) $`\mathrm{\Delta }_0<0`$ such that at $`T=0K`$ all the available electrons form local pairs (of ”d” electrons) (the ”local pair” regime, $`2n_Bn_c`$) (LP);
(ii) $`\mathrm{\Delta }_0>0`$ such that the ”c” electron band is filled up to the Fermi level $`\mu =\mathrm{\Delta }_0`$ and the remaining electrons are in the form of local pairs (the c+d regime or Mixed, $`0<2n_B,n_c<2)`$ (LP+E);
(iii) $`\mathrm{\Delta }_0>0`$ such that the Fermi level $`\mu <\mathrm{\Delta }_0`$ and consequently at $`T=0`$K all the available electrons occupy the ”c” electron states (the c-regime or”BCS”, $`n_c2n_B`$) (E).

For $`|I_0|0`$, in the case (ii) superconductivity is due to the interchange between local pairs of ”d” electrons and pairs of ”c” electrons. In this process ”c” electrons become ”polarized” into Cooper pairs and ”d” electron pairs increase their mobility by decaying into ”c” electron pairs. To this intermediate case neither the standard BCS picture nor the picture of local pairs applies and superconductivity has a ”mixed” character. The system shows features which are intermediate between the BCS and preformed local pair regime. This regards the energy gap in the single-electron excitation spectrum ($`E_g(0))`$, the $`k_BT_c/E_g(0)`$ ratio, the critical fields, the Ginzburg ratio $`\kappa `$, the width of the critical regime as well as the normal state properties. In case (i) the local pairs of ”d” electrons can move via a mechanism of virtual excitations into empty c-electrons states. Such a mechanism gives rise to the long range hopping of pairs of ”d” electrons (in analogy to the RKKY interaction for s-d mechanism in the magnetic equivalent). The superconducting properties are analogous to those of a pure local pair superconductor . In case (iii), on the contrary, we find a situation which is similar to the BCS case: pairs of ”c” electrons with opposite momenta and spins are exchanged via virtual transitions into local pair states.

The generic phase diagrams for different pairing symmetries plotted as a function of the position of the LP level $`\mathrm{\Delta }_0`$ at fixed $`n`$ are shown in Fig.2a-c. In all the cases one observes a drop in the superfluid stiffness (and in the KT transition temperature) when the bosonic level reaches the bottom of the $`c`$-electron band and the system approaches the LP limit. In the opposite, BCS like limit, $`T_c^{KT}`$ approaches asymptotically $`T_c^{MFA}`$, with a narrow fluctuation regime. Between the KT and MFA temperatures, the phase fluctuation effects are important. In this regime a pseudogap in the c-electron spectrum will develop and the normal state of LP and itinerant fermions can exhibit non-Fermi liquid properties .

The analysis of the Mixed-LP crossover indicates that when the LP level is lowered and reaches the bottom of the fermionic band an effective attraction between fermions becomes strong, since it varies as $`I^2/(2\mathrm{\Delta }_02\mu )`$ and $`\mu \mathrm{\Delta }_0`$ (Eq. 10). In this regime the density of $`c`$ electrons is low and formation of bound $`c`$-electron pairs occurs. It gives rise to an energy gap in the single-electron spectrum independently of the pairing symmetry. We evaluated the binding energies of $`c`$\- electron pairs and found that $`T_c^{MFA}`$ scales with the half of their binding energy for $`\mathrm{\Delta }_00`$. The superconducting transition temperature is here always much lower than the $`c`$-pair formation temperature ($`T_c^{MFA}`$) and rapidly decreases with $`|\mathrm{\Delta }_0/D|`$. In such a case, the superconducting state can be formed by two types of coexisting (composite) bosons: preformed $`c`$\- electron pairs and LP’s .

Comparing $`T_c`$ vs $`\mathrm{\Delta }_0`$ plots for various pairings one finds that in the case of nn hopping only, the $`d`$ and $`s`$ -wave pairings are favorable for higher concentration of $`c`$-electrons, while the $`s^{}`$-wave can be stable at low $`n_c`$. This is clearly seen from Fig.3, where the plots of $`T_c^{MFA}`$ vs $`\mathrm{\Delta }_0`$ for $`s,s^{}`$ and $`d`$-pairings are given for $`n=1,|I_0|/D=0.25`$, and from Fig.4 which shows variation of $`n_c`$ and $`n_B`$ with $`\mathrm{\Delta }_0/D`$ for the same values of $`n`$ and $`\mathrm{\Delta }_0/D`$ for $`d`$ -wave pairing (for $`s`$ and $`s^{}`$ -wave symmetries the plots of $`n_c`$ and $`n_B`$ are very similar).

A strongly nonmonotonous behavior of $`T_c`$ vs $`\mathrm{\Delta }_0`$ for $`s^{}`$-wave pairing and $`n>1`$ (cf. Fig.2b) reflects the fact that this type of pairing in a single band model is realized for low electron (or hole densities), if $`t_2=0`$ . In general, with increasing $`n_c`$ $`(\mathrm{\Delta }_0/D`$) one can have a sequence of transitions: $`s^{}d`$ (or $`s`$)$`s^{}`$ (if $`1<n2`$), or the transitions: $`s^{}d`$ (or $`s`$) (if $`n1`$).
The nnn hopping $`t_2`$ (with opposite sign to $`t`$) can strongly enhance $`T_c`$ for $`d`$-wave symmetry, moreover it favors the $`d`$ and $`s`$-wave pairings for lower values of $`n_c`$ .

The region between $`T_c^{MFA}`$ and $`T_c^{KT}`$, where the system can exhibit a pseudogap, expands with increasing intersubsystem coupling $`|I_0|`$. It is seen from Fig.5 that, except for $`|I_0|/D1`$ the coupling dependences of $`T_c^{MFA}`$ and $`T_c^{KT}`$ are qualitatively different. $`T_c^{MFA}`$ is an increasing function of $`|I_0|`$ for all the pairing symmetries. On the other hand, $`T_c^{KT}`$ vs $`|I_0|`$ increases first, goes through a round maximum and then decreases (similarly as it is observed in the attractive Hubbard model). The position of the maximum corresponds to the intermediate values of $`|I_0|/D`$ and it depends on the pairing symmetry as well as the values of $`\mathrm{\Delta }_0/D`$ and $`n`$. For large $`|I_0|`$, the $`T_c^{KT}`$ are close to the upper bound for the phase ordering temperature which is given by $`\pi \rho _s(0)/2`$.

As for the evolution of superconducting properties with increasing $`n`$ one finds three possible types of density driven changeovers : (i) for $`2\mathrm{\Delta }_0/D0`$, ”BCS”$``$ Mixed $``$ ”BCS”; (ii) for $`\mathrm{\Delta }_0/D>2`$: ”BCS” $``$ ”LP” and (iii) for $`\mathrm{\Delta }_0/D<0`$: ”LP” $``$ ”BCS”. Only if the LP level is deeply located below the bottom of $`c`$-band, the system remains in the LP regime for any $`n2`$.

Some of our findings can be qualitatively related to experimental results for the cuprate HTS where a pseudogap exists. It has been suggested by ARPES experiments, that for underdoped cuprates the Fermi surface in the pseudogap phase is truncated around the corners due to the formation of preformed (bosonic) pairs with charge $`2e`$, whereas the ”electrons” on the diagonals remain unpaired \[7a\]. In the present two component model such a situation is obtained when LP’s and $`c`$-electrons coexist in the mixed regime. The linear $`T`$-dependence of the superfluid density has been observed experimentally in copper oxides and also in several organic superconductors. This points to an order parameter of $`d_{x^2y^2}`$-wave symmetry and existence of nodal quasiparticles. In the present model the gap ratio is nonuniversal and can deviate strongly from BCS predictions (particularly in the d-wave case for which it is always enhanced) . This feature is also found for several exotic superconductors. The Uemura plots i.e., the $`T_c`$ vs zero-temperature phase stiffness $`\rho _s(0)`$ (and the scaling $`T_c\rho _s(0)`$) reported for cuprates and organic superconductors can be reproduced within the model for extended $`s`$\- and $`d`$-wave order parameter symmetry (, Figs.11,13 in Ref.).
Finally, we point out that to relate more quantitatively the present two-component scenario to the phase diagram of cuprates one has to consider doping dependences of the model parameters.

## Acknowledgements

This work was supported in part by the Polish State Committee for Scientific Research (KBN), Project No. 2 P03B 154 22.