1 Introduction

In the studies of a variety of the strong interaction processes a key role is played by the hadronic vacuum polarization function \(\varPi (q^2)\), the related function R(s), and the Adler function \(D(Q^2)\). In particular, these functions govern such processes as the electron–positron annihilation into hadrons, inclusive hadronic decays of \(\tau \) lepton and Z boson, as well as the hadronic contributions to precise electroweak observables, such as the muon anomalous magnetic moment \((g-2)_{\mu }\) and the running of the electromagnetic fine structure constant. The theoretical analysis of these processes constitutes a decisive self-consistency test of Quantum Chromodynamics (QCD) and entire Standard Model, that, in turn, puts robust restrictions on a possible New Physics beyond the latter. Additionally, the energy scales relevant to the foregoing strong interaction processes span from the infrared to ultraviolet domain, so that their theoretical investigation provides a native framework for a profound study of both perturbative and intrinsically nonperturbative aspects of hadron dynamics. It is worth noting also that a majority of the aforementioned processes are of a direct relevance to the physics at the currently designed Future Collider Projects, such as the Future Circular Collider FCC-ee [1], Circular Electron–Positron Collider CEPC (its first phase) [2], the International Linear Collider ILC [3], the Compact Linear Collider CLIC [4], as well as the E989 experiment at Fermilab [5], the E34 experiment at J-PARC [6], and others.

In fact, over the past decades the perturbative approach to QCD remains a basic tool for the theoretical exploration of the hadronic physics. However, the QCD perturbation theory can be directly applied to the study of the strong interaction processes only in the spacelike (Euclidean) domain, whereas the proper description of hadron dynamics in the timelike (Minkowskian) domain additionally requires the pertinent dispersion relations. Specifically, the dispersion relation for the R-ratio of electron–positron annihilation into hadrons converts the physical kinematic restrictions on the process on hand into the mathematical form and determines the way how the “timelike” observable R(s) is related to the “spacelike” quantity \(D(Q^2)\), the corresponding perturbative input being embodied by the so-called spectral function. Since the calculation of the latter at the higher-loop levels constitutes a rather challenging task, one commonly resorts to an approximate form of the R-ratio, namely, its truncated re-expansion at high energies. At the same time, one has to be aware that at any given loop level such re-expansion generates an infinite number of the so-called \(\pi ^2\)-terms (which may not necessarily be small enough to be safely discarded at the higher orders), which also worsen the loop convergence of the resulting approximate R-ratio.

The primary objective of the paper is to explicitly derive a general form of the spectral function required for the proper evaluation of the R-ratio and to study its properties up to the five-loop level. It is also of an apparent interest to calculate the R-ratio itself, to examine its higher-loop convergence, and to elucidate the impact of the omitted higher-order \(\pi ^2\)-terms on its truncated re-expanded approximation.

The layout of the paper is as follows. In Sect. 2 the essentials of continuation of the spacelike perturbative results into the timelike domain are expounded. In Sect. 3 the one-loop expression for the R-ratio is explicated and its approximations are discussed. In Sect. 4 the derivation of a general form of the commonly employed approximate expression for the R-ratio (which constitutes its truncated re-expansion at high energies) is delineated, the appearance of the pertinent \(\pi ^2\)-terms is elucidated, and their basic features are studied. In Sect. 5 the explicit form of the spectral function required for the proper calculation of the R-ratio is obtained and its properties at the higher-loop levels are examined. By making use of the derived spectral function in Sect. 6 the proper expression for the R-ratio is calculated up to the five-loop level and its properties are studied. Additionally, the obtained R-ratio is juxtaposed with its commonly employed approximation and the impact of the omitted higher-order \(\pi ^2\)-terms on the latter is discussed. In Sect. 7 the basic results are summarized.

2 R-ratio of electron–positron annihilation into hadrons

As noted above, the theoretical analysis of certain strong interaction processes relies on the hadronic vacuum polarization function \(\varPi (q^2)\), which is defined as the scalar part of the hadronic vacuum polarization tensor,

$$\begin{aligned} \varPi _{\mu \nu }(q^2)= & {} i \int \mathrm{d}^4x\,e^{i q x} \bigl \langle 0 \bigl |\, T \left\{ J_{\mu }(x)\, J_{\nu }(0)\right\} \bigr | 0 \bigr \rangle \nonumber \\= & {} \frac{i}{12\pi ^2} (q_{\mu }q_{\nu } - g_{\mu \nu }q^2) \varPi (q^2). \end{aligned}$$
(1)

For the processes involving final state hadrons the function \(\varPi (q^2)\) (1) has the only cut along the positive semiaxis of real \(q^2\) starting at the hadronic production threshold \(q^2 \ge 4m_{\pi }^2\) (the discussion of this issue can be found in, e.g., Ref. [7]). In particular, the Feynman amplitude of the respective process vanishes for the energies below the threshold, which expresses the physical fact that the production of the final state hadrons is kinematically forbidden for \(q^2 < 4m_{\pi }^2\). In turn, the known location of the cut of function \(\varPi (q^2)\) in the complex \(q^2\) plane enables one to write down the pertinent dispersion relation

$$\begin{aligned} \varDelta \varPi (q^2,\, q_0^2) = (q^2 - q_0^2) \int \limits _{4m_{\pi }^2}^{\infty } \frac{R(\sigma )}{(\sigma -q^2)(\sigma -q_0^2)}\, \mathrm{d}\sigma , \end{aligned}$$
(2)

with the once-subtracted Cauchy integral formula being employed. In Eq. (2) \(\varDelta \varPi (q^2,\, q_0^2) = \varPi (q^2) - \varPi (q_0^2)\), whereas R(s) stands for the discontinuity of the hadronic vacuum polarization function across the physical cut

$$\begin{aligned} R(s) = \frac{1}{2 \pi i} \lim _{\varepsilon \rightarrow 0_{+}} \varDelta \varPi (s+i\varepsilon ,s-i\varepsilon ). \end{aligned}$$
(3)

This function is commonly identified with the so-called R-ratio of electron–positron annihilation into hadrons \(R(s) = \sigma (e^{+}e^{-} \rightarrow \text {hadrons}; s)/\sigma (e^{+}e^{-} \rightarrow \mu ^{+}\mu ^{-}; s)\), with \(s=q^2>0\) being the timelike kinematic variable, namely, the center-of-mass energy squared.

In practice one deals with the Adler function \(D(Q^2)\) [8], which is defined as the logarithmic derivative of the hadronic vacuum polarization function (1)

$$\begin{aligned} D(Q^2) = - \frac{\mathrm{d}\, \varPi (-Q^2)}{\mathrm{d} \ln Q^2}, \end{aligned}$$
(4)

with \(Q^2 = -q^2 > 0\) being the spacelike kinematic variable. Note that the subtraction point \(q_0^2\) entering Eq. (2) does not appear in Eqs. (3) and (4). The widely employed dispersion relation for the Adler function follows immediately from Eqs. (2) and (4), specifically [8],

$$\begin{aligned} D(Q^2) = Q^2 \int \limits _{4m_{\pi }^2}^{\infty } \frac{R(\sigma )}{(\sigma +Q^2)^2}\,\mathrm{d}\sigma . \end{aligned}$$
(5)

In particular, this dispersion relation enables one to extract the experimental prediction for the Adler function by making use of the corresponding experimental data on electron–positron annihilation into hadrons. However, to obtain the theoretical expression for the R-ratio itself, the relation inverse to Eq. (5) is required. The latter can be obtained by integrating Eq. (4) in finite limits, which yields [9,10,11]

$$\begin{aligned} R(s) = \frac{1}{2 \pi i} \lim _{\varepsilon \rightarrow 0_{+}} \int \limits _{s + i \varepsilon }^{s - i \varepsilon } D(-\zeta )\,\frac{\mathrm{d} \zeta }{\zeta }. \end{aligned}$$
(6)

Specifically, this equation relates the R-ratio to the theoretically calculable Adler function and provides a native way to properly account for the effects due to continuation of the spacelike perturbative results into the timelike domain. The integration contour in Eq. (6) lies in the region of analyticity of the integrand; see Fig. 1. Note also that the relation which expresses the hadronic vacuum polarization function in terms of the Adler function can be obtained in a similar way, namely [12]

$$\begin{aligned} \varDelta \varPi (-Q^2,\, -Q_0^2) = - \int \limits _{Q_0^2}^{Q^2} D(\zeta ) \frac{\mathrm{d} \zeta }{\zeta }, \end{aligned}$$
(7)

where \(Q^{2}\) and \(Q_{0}^{2}\) stand for the spacelike kinematic variable and the subtraction point, respectively. Basically, Eqs. (2)–(7) constitute the complete set of relations, which express the functions on hand in terms of each other. It is worth mentioning here that, in general, the pattern of applications of the dispersion relations in theoretical particle physics is quite diverse. For example, among the latter are such issues as the refinement of chiral perturbation theory [13,14,15,16,17,18,19,20], the accurate determination of parameters of resonances [21,22,23,24], the assessment of the hadronic light-by-light scattering [25,26,27,28,29,30], as well as many others.

Fig. 1
figure 1

The integration contour in Eq. (6). The physical cut \(\zeta \ge 4m_{\pi }^2\) of the Adler function \(D(-\zeta )\) (4) is shown along the positive semiaxis of real \(\zeta \)

It is necessary to outline that the derivation of dispersion relations (2)–(7) is based only on the kinematics of the process on hand and involves neither model-dependent phenomenological assumptions nor additional approximations. In turn, the relations (2)–(7) impose a number of strict physical inherently nonperturbative constraints on the functions \(\varPi (q^2)\), R(s), and \(D(Q^2)\), which should definitely be taken into account when one intends to go beyond the limits of applicability of the QCD perturbation theory. It is worthwhile to note that these nonperturbative restrictions have been merged with the corresponding perturbative input in the framework of dispersively improved perturbation theory (DPT) [31,32,33,34] (its preliminary formulation was discussed in Refs. [35,36,37,38]). In particular, the DPT enables one to overcome some intrinsic difficulties of the QCD perturbation theory and to extend its applicability range towards the infrared domain; see Ref. [31] and the references therein for the details.

In the framework of perturbation theory the Adler function (4) takes the form of the power series in the so-called QCD couplant \(a^{(\ell )}_{\text {s}}(Q^2) = \alpha ^{(\ell )}_{\text {s}}(Q^2)\, \beta _{0}/(4\pi )\), namely

$$\begin{aligned}&D^{(\ell )}_{\mathrm{pert}}(Q^2) = 1 + d^{(\ell )}_{\mathrm{pert}}(Q^2),\nonumber \\&\quad d^{(\ell )}_{\mathrm{pert}}(Q^2) = \sum _{j=1}^{\ell } d_{j} \left[ a^{(\ell )}_{\text {s}}(Q^2) \right] ^{j}. \end{aligned}$$
(8)

In this equation \(\ell \) specifies the loop level, \(d_1 = 4/\beta _0\), \(\beta _0 = 11 - 2n_{f}/3\), \(n_{f}\) is the number of active flavors, and the common prefactor \(N_{\text {c}}\sum _{f=1}^{n_{f}} Q_{f}^{2}\) is omitted throughout, where \(N_{\text {c}}=3\) denotes the number of colors and \(Q_{f}\) stands for the electric charge of fth quark. The QCD couplant \(a^{(\ell )}_{\text {s}}(Q^2)\) entering Eq. (8) can be represented as the double sum

$$\begin{aligned} a^{(\ell )}_{\text {s}}(Q^2) = \sum _{n=1}^{\ell }\sum _{m=0}^{n-1} b^{m}_{n}\, \frac{\ln ^{m}(\ln z)}{\ln ^n z}, \end{aligned}$$
(9)

where \(z=Q^2/\varLambda ^2\) and \(b^{m}_{n}\) (the integer superscript m is not to be confused with respective power) stands for the combination of the \(\beta \) function perturbative expansion coefficients, specifically, \(b^{0}_{1}=1\), \(b^{0}_{2}=0\), \(b^{1}_{2}=-\beta _{1}/\beta ^{2}_{0}\), etc. The Adler function perturbative expansion coefficients \(d_{j}\) were calculated up to the four-loop level (\(1 \le j \le 4\)); see Refs. [39,40,41] and the references therein, whereas for the five-loop coefficient \(d_{5}\) only numerical estimation [42] is available so far. The numerical values of the perturbative coefficients \(d_{j}\) (8) are listed in Table 1. In turn, the \(\beta \) function perturbative expansion coefficients \(\beta _{j}\) have been calculated up to the five-loop level (\(0 \le j \le 4\)); see Refs. [43, 44] and the references therein for the details.

Table 1 Numerical values of the Adler function perturbative expansion coefficients \(d_j\) (8). In the last column the numerical estimation of the five-loop coefficient \(d_5\) [42] is listed

In what follows the nonperturbative aspects of the strong interactions will be disregarded and a primary attention will be given to the theoretical description of the R-ratio of electron–positron annihilation into hadrons at moderate and high energies. For this purpose the effects due to the masses of the involved particles can be safely neglected (a discussion of the impact of such effectsFootnote 1 on the low-energy behavior of the functions \(\varPi (q^2)\), R(s), and \(D(Q^2)\) can be found in, e.g., Refs. [31,32,33,34, 45,46,47,48,49,50,51,52,53]). Additionally, for the scheme-dependent perturbative coefficients \(\beta _{j}\) and \(d_{j}\) the \(\overline{\mathrm{MS}}\)-scheme will be assumed and for the uncalculated yet five-loop coefficient \(d_5\) its numerical estimation [42] will be employed.

Thus, in the massless limit the relation (6) can be represented as (see also Refs. [54, 55])

$$\begin{aligned} R^{(\ell )}(s) = 1 + r^{(\ell )}(s), \quad r^{(\ell )}(s) = \int \limits _{s}^{\infty } \rho ^{(\ell )}(\sigma )\, \frac{\mathrm{d} \sigma }{\sigma }, \end{aligned}$$
(10)

where

$$\begin{aligned} \rho ^{(\ell )}(\sigma ) = \frac{1}{2 \pi i} \lim _{\varepsilon \rightarrow 0_{+}} \Bigl [d^{(\ell )}(-\sigma - i \varepsilon ) - d^{(\ell )}(-\sigma + i \varepsilon ) \Bigr ]\nonumber \\ \end{aligned}$$
(11)

stands for the spectral function and \(d^{(\ell )}(Q^2)\) denotes the \(\ell \)-loop strong correction to the Adler function. As mentioned above, only perturbative contributions will be retained in Eq. (11) hereinafter, which makes Eq. (10) identical to that of both the foregoing DPT [31,32,33,34] and the so-called analytic perturbation theoryFootnote 2 (APT) [54,55,56,57]. It has to be noted that, in general, the perturbative spectral function at small values of its argument may be altered by the terms of an intrinsically nonperturbative nature. For instance, the nonperturbative models discussed in Refs. [105,106,107,108,109,110,111,112] constitute a superposition of the perturbative spectral function with the so-called “flat” terms (which by definition do not affect the corresponding perturbative results at high energies), whereas the models [113, 114] modify the low-energy behavior of the perturbative spectral function proceeding from certain phenomenological assumptions.

It is worthwhile to mention also that a “naive” approach to continue the spacelike perturbative result (8) into the timelike domain consists in merely identifying the timelike kinematic variable (\(s=q^2\)) with the spacelike one (\(Q^2=-q^2\)), i.e.,

$$\begin{aligned} R^{(\ell )}_{\text {naive}}(s) = D^{(\ell )}_{\text {pert}}(|s|) = 1 + \sum _{j=1}^{\ell } d_{j} \left[ a^{(\ell )}_{\text {s}}(|s|) \right] ^{j}. \end{aligned}$$
(12)

However, as thoroughly discussed in Refs. [115,116,117], this prescription yields a misleading result, which differs from the proper one (10) even in the deep ultraviolet asymptotic, see also Refs. [42, 118, 119] as well as [31] and the references therein.

3 R-ratio at the one-loop level

Let us address now the R-ratio of electron–positron annihilation into hadrons at the one-loop level. As discussed in the previous section, for this purpose the spectral function \(\rho (\sigma )\) (11), which enters the pertinent integral representation (10), is required. Since the involved strong correction to the Adler function (8) takes a simple form at the one-loop level,

$$\begin{aligned} d^{(1)}_{\text {pert}}(Q^2) = d_{1}\,a^{(1)}_{\text {s}}(Q^2), \quad a^{(1)}_{\text {s}}(Q^2) = \frac{1}{\ln (Q^2/\varLambda ^2)}, \end{aligned}$$
(13)

and

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0_{+}} \ln (x \pm i\varepsilon ) = \ln |x| \pm i\pi \theta (-x), \end{aligned}$$
(14)

the calculation of \(\rho ^{(1)}(\sigma )\) (11) appears to be quite straightforward. Specifically, the one-loop spectral function (11) for the positive values of its argument reads

$$\begin{aligned} \rho ^{(1)}(\sigma ) = d_{1}\,\frac{1}{2 \pi i} \left[ \frac{1}{\ln (\sigma /\varLambda ^2) - i \pi } - \frac{1}{\ln (\sigma /\varLambda ^2) + i \pi }\right] ,\nonumber \\ \end{aligned}$$
(15)

which, in turn, can be represented as

$$\begin{aligned} \rho ^{(1)}(\sigma ) = d_{1} {\bar{\rho }}^{(1)}_{1}(\sigma ), \; {\bar{\rho }}^{(1)}_{1}(\sigma ) = \frac{1}{y^2+\pi ^2}, \; y=\ln \left( \frac{\sigma }{\varLambda ^2} \right) .\nonumber \\ \end{aligned}$$
(16)

In these equations \(\theta (x)\) is the Heaviside unit step function [i.e., \(\theta (x)=1\) if \(x \ge 0\) and \(\theta (x)=0\) otherwise], \(d_{1}=4/\beta _{0}\), and \(\beta _{0}=11 - 2n_{f}/3\). The plot of the one-loop spectral function \({\bar{\rho }}^{(1)}_{1}(\sigma )\) (16) is displayed in Fig. 2. As one can infer from this figure, the function on hand assumes the values in the interval \(0 \le {\bar{\rho }}^{(1)}_{1}(\sigma ) \le 1/\pi ^2\) and decreases as \(1/y^2\) in both ultraviolet (\(y\rightarrow \infty \)) and infrared (\(y\rightarrow -\infty \)) asymptotics.

Fig. 2
figure 2

The one-loop spectral function \({\bar{\rho }}^{(1)}_{1}(\sigma )\) (16)

Fig. 3
figure 3

Plot A: The one-loop “timelike” effective couplant \(a^{(1)}_{\text { TL}}(s)\) [Eq. (17), solid curve] and the naive continuation of the one-loop perturbative couplant into the timelike domain \(a^{(1)}_{\text {s}}(|s|)\) [Eq. (13), dashed curve]. Plot B: The one-loop “timelike” effective couplant \(a^{(1)}_{\text { TL}}(s)\) [Eq. (17), solid curve] and its re-expansions at various energy intervals: \(\ln w < -\pi \) [Eq. (18), dotted curves], \(-\pi< \ln w < \pi \) [Eq. (19), dot-dashed curves], and \(\ln w > \pi \) [Eq. (20), dashed curves]. Numerical labels indicate the highest absolute value of the power of \(\ln w\) retained in the re-expansions on hand. The boundaries of the convergence ranges of Eqs. (18)–(20) are marked by vertical solid lines

Then the corresponding one-loop strong correction to the R-ratio can also easily be obtained in an explicit form. Specifically, the integration (10) of the one-loop spectral function (16) yieldsFootnote 3

$$\begin{aligned} r^{(1)}(s) = d_{1}a^{(1)}_{\text { TL}}(s), \quad a^{(1)}_{\text { TL}}(s)= \frac{1}{2} - \frac{1}{\pi }\arctan \left( \frac{\ln w}{\pi }\right) ,\nonumber \\ \end{aligned}$$
(17)

where \(w=s/\varLambda ^2\). The function \(a^{(1)}_{\text { TL}}(s)\) (17) constitutes the one-loop couplant, which properly accounts for the effects due to continuation of the spacelike perturbative expression (13) into the timelike domain. It is worthwhile to note here that Eq. (17) has first appeared in Ref. [120] and only afterwards was derived in Refs. [9, 10, 12, 54, 55].

Figure 3A displays the one-loop “timelike” effective couplant \(a^{(1)}_{\text { TL}}(s)\) (17) and the “naive” continuation of the one-loop perturbative couplant \(a^{(1)}_{\text {s}}(Q^2)\) (13) into the timelike domain (12). As one can infer from this figure, at high  energies the two couplants approach each other. At the same time, at moderate energies the deviation between the functions on hand becomes significant, whereas in the infrared domain their behavior turns out to be qualitatively different. Specifically, the function \(a^{(1)}_{\text {s}}(|s|)\) (13) diverges at low energies due to the infrared unphysical singularities, whereas \(a^{(1)}_{\text { TL}}(s)\) (17) is a smooth monotone decreasing function of its argument, which contains no singularities for \(s{>}0\).

A somewhat simpler but approximate form of the strong correction to the R-ratio can be obtained by its further re-expansion. Specifically, for this purpose one splits the entire energy range \(0< s < \infty \) into three intervals (namely, \(\ln w < -\pi \), \(-\pi< \ln w < \pi \), and \(\ln w> \pi \)) and applies the Taylor expansion to r(s) in each of those intervals. At the one-loop level the implementation of these steps for the expression (17) yields

$$\begin{aligned} a^{(1)}_{\text { TL}}(s)\simeq & {} 1 + \frac{1}{\ln w} - \frac{1}{3}\frac{\pi ^2}{\ln ^{3}w} +\mathcal {O}\left( \frac{1}{\ln ^{5}w} \right) ,\nonumber \\&\quad \ln w < -\pi , \end{aligned}$$
(18)
$$\begin{aligned} a^{(1)}_{\text { TL}}(s)\simeq & {} \frac{1}{2} - \frac{\ln w}{\pi ^2} + \frac{1}{3}\frac{\ln ^{3}w}{\pi ^4} +\mathcal {O}( \ln ^{5}w ),\nonumber \\&\quad -\pi< \ln w < \pi , \end{aligned}$$
(19)
$$\begin{aligned} a^{(1)}_{\text { TL}}(s)\simeq & {} \frac{1}{\ln w} - \frac{1}{3}\frac{\pi ^2}{\ln ^{3}w} +\mathcal {O}\left( \frac{1}{\ln ^{5}w} \right) , \quad \ln w > \pi ,\nonumber \\ \end{aligned}$$
(20)

where \(w=s/\varLambda ^2\). In particular, as one can infer from Fig. 3B, the re-expansions (18)–(20) may provide an accurate approximation of the function (17) in the aforementioned energy intervals, if the number of retained terms is large enough. However, as one can note, the convergence of the re-expansions (18)–(20) becomes worse when the energy scale approaches the delimiting values \(\sqrt{s}/\varLambda = \exp (\pm \pi /2)\); see also discussion of this issue in Sect. 6.

In fact, there is another equivalent way to obtain the re-expansion of the strong correction to the R-ratio at high energies. Specifically, instead of following the lines described above, one can expand the corresponding spectral function \(\rho ^{(\ell )}(\sigma )\) (11) and then perform the integration in Eq. (10). In particular, at the one-loop level the Taylor expansion of \({\bar{\rho }}^{(1)}_{1}(\sigma )\) (16) for \(\sqrt{\sigma }/\varLambda >\exp (\pi /2)\) reads

$$\begin{aligned} {\bar{\rho }}^{(1)}_{1}(\sigma ) \simeq \frac{1}{y^2} - \frac{\pi ^2}{y^4} + \mathcal {O}\left( \frac{1}{y^6}\right) ,\quad y=\ln \left( \frac{\sigma }{\varLambda ^2}\right) , \end{aligned}$$
(21)

which, after its integration in Eq. (10), yields the result identical to Eq. (20). It is the latter prescription that will be employed in the next section for the derivation of an approximate expression for the R-ratio at high energies at an arbitrary loop level.

4 Re-expansion of the R-ratio at high energies: \(\pi ^2\)-terms

As outlined in Sect. 2, the strong correction to the R-ratio of electron–positron annihilation into hadrons (10) can be represented as

$$\begin{aligned} r(s) = \int \limits _{s}^{\infty } \rho (\sigma )\,\frac{\mathrm{d} \sigma }{\sigma } = \int \limits _{\ln w}^{\infty }\rho _{y}(y)\, \mathrm{d}y, \quad w=\frac{s}{\varLambda ^2}, \end{aligned}$$
(22)

where \(y=\ln (\sigma /\varLambda ^2)\) and \(\rho _{y}(y) = \rho [\varLambda ^2\exp (y)]\) denotes the corresponding spectral function (11)

$$\begin{aligned} \rho _{y}(y) = \frac{1}{2 \pi i}\, \Bigl [d_{y}(y - i \pi ) - d_{y}(y + i \pi )\Bigr ], \end{aligned}$$
(23)

with Eq. (14) being employed. In Eq. (23) \(d_{y}(y) = d\bigl [\varLambda ^2\exp (y)\bigr ]\) stands for the strong correction to the Adler function being expressed in terms of \(y=\ln (\sigma /\varLambda ^2)\). Applying to the latter the Taylor expansion

$$\begin{aligned} d_{y}(y \pm i \pi ) = d_{y}(y) + \sum _{n=1}^{\infty } \frac{(\pm i\pi )^n}{n !}\frac{\mathrm{d}^{n}}{\mathrm{d} y^{n}}\, d_{y}(y), \quad |y| > \pi ,\nonumber \\ \end{aligned}$$
(24)

one can approximateFootnote 4 the spectral function (23) by

$$\begin{aligned} \rho _{y}(y) = - \frac{\mathrm{d}}{\mathrm{d}y}\, d_{y}(y) - \sum _{n=1}^{\infty } \frac{(-1)^{n} \pi ^{2n}}{(2n+1)!}\, \frac{\mathrm{d}^{2n+1}}{\mathrm{d}y^{2n+1}}\, d_{y}(y). \end{aligned}$$
(25)

Therefore, for \(\sqrt{s}/\varLambda > \exp (\pi /2) \simeq 4.81\) the strong correction to the R-ratio (22) acquires the following form:

$$\begin{aligned} r(s) = d(|s|) + \sum _{n=1}^{\infty } \frac{(-1)^{n} \pi ^{2n}}{(2n+1)!}\, \frac{\mathrm{d}^{2n}}{\mathrm{d} y^{2n}}\, d_{y}(y) \Biggr |_{y=\ln w}. \end{aligned}$$
(26)

In particular, this equation implies that the strong correction to the R-ratio, being re-expanded at high energies, reproduces the naive continuation of the Adler function into the timelike domain [the first term on the right-hand side of Eq. (26)] and additionally produces an infinite number of the so-called \(\pi ^{2}\)-terms.

Then at the \(\ell \)-loop level the perturbative expression for the strong correction to the Adler function reads (8)

$$\begin{aligned} d^{(\ell )}_{y}(y) = \sum _{j=1}^{\ell } d_{j} \left[ a^{(\ell )}_{y}(y)\right] ^{j}, \end{aligned}$$
(27)

where \(a^{(\ell )}_{y}(y)=a^{(\ell )}_{\text {s}}\bigl [\varLambda ^2\exp (y)\bigr ]\) is the \(\ell \)-loop perturbative couplant being expressed in terms of \(y=\ln (\sigma /\varLambda ^2)\). Since the latter satisfies the renormalization group equation

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}y}\, a^{(\ell )}_{y}(y) = - \sum _{j=0}^{\ell -1} B_j \left[ a^{(\ell )}_{y}(y)\right] ^{j+2}, \quad B_j=\frac{\beta _j}{\beta ^{j+1}_{0}}, \end{aligned}$$
(28)

the nth derivative of the jth power of the \(\ell \)-loop couplant takes the form

$$\begin{aligned}&\frac{\mathrm{d}^{n}}{\mathrm{d}y^{n}} \left[ a^{(\ell )}_{y}(y)\right] ^{j} =(-1)^{n} \sum _{k_{1}=0}^{\ell -1} \ldots \sum _{k_{n}=0}^{\ell -1} \left[ a^{(\ell )}_{y}(y) \right] ^{j+n+k_{1}+\cdots +k_{n}}\nonumber \\&\quad \times \left( \prod _{p=1}^{n}B_{k_{p}} \right) \left[ \prod _{t=0}^{n-1}(j+t+k_{1}+\cdots +k_{t})\right] . \end{aligned}$$
(29)

Thus, at high energies the \(\ell \)-loop strong correction to the R-ratio (22) can be approximated by (see also Ref. [121])

$$\begin{aligned}&r^{(\ell )}(s)=\sum _{j=1}^{\ell } d_{j} \left[ a^{(\ell )}_{\mathrm{s}}(|s|)\right] ^{j} - \sum _{j=1}^{\ell } d_{j} \sum _{n=1}^{\infty } \frac{(-1)^{n+1}\pi ^{2n}}{(2n+1)!}\nonumber \\&\quad \times \sum _{k_{1}=0}^{\ell -1} \cdots \sum _{k_{2n}=0}^{\ell -1} \left( \prod _{p=1}^{2n}B_{k_{p}} \right) \nonumber \\&\quad \times \left[ \prod _{t=0}^{2n-1} (j+t+k_{1}+k_{2}+\cdots +k_{t})\right] \nonumber \\&\quad \times [a^{(\ell )}_{\mathrm{s}}(|s|)]^{j+2n+k_{1}+k_{2}+\cdots +k_{2n}},\quad \frac{\sqrt{s}}{\varLambda } > \exp \left( \frac{\pi }{2}\right) . \end{aligned}$$
(30)

The obtained re-expansion of the strong correction to the R-ratio (30) constitutes the sum of naive continuation of the strong correction to the Adler function into the timelike domain (12) and an infinite number of the \(\pi ^2\)-terms. Equation (30) explicitly proves the fact that at any given loop level the re-expansion of the strong correction to the R-ratio at high energies can be reduced to the form of power series in the naive continuation of the perturbative couplant into the timelike domain \(a^{(\ell )}_{\text {s}}(|s|)\). As one can also note, in the re-expansion (30) the coefficients \(d_j\) corresponding to various orders of perturbation theory turn out to be all mixed up, i.e., the \(\ell \)-loop contribution to Eq. (10) appears to be re-distributed over the higher-order terms.

As discussed earlier, if the number of terms retained in Eq. (30) is large enough, then it can provide a rather accurate approximation of the strong correction to the R-ratio (10) for \(\sqrt{s}/\varLambda > \exp (\pi /2) \simeq 4.81\). However, one usually truncates the re-expansion (30) at the order \(\ell \), thereby neglecting all the higher-order \(\pi ^2\)-terms,Footnote 5 which results in the following expression commonly employed in practical applications:

$$\begin{aligned} R^{(\ell )}_{\mathrm{appr}}(s) = 1 + r^{(\ell )}_{\mathrm{appr}}(s), \quad r^{(\ell )}_{\mathrm{appr}}(s) = \sum _{j=1}^{\ell } r_{j} \left[ a^{(\ell )}_{\mathrm{s}}(|s|)\right] ^{j},\nonumber \\ \end{aligned}$$
(31)

where

$$\begin{aligned} r_{j} = d_{j} - \delta _{j}. \end{aligned}$$
(32)

In this equation the \(d_j\) stand for the Adler function perturbative expansion coefficients (8), whereas the \(\delta _{j}\) embody the contributions of the relevant \(\pi ^2\)-terms.

Equation (30) implies that the \(\pi ^2\)-terms do not appear in the first and second orders of perturbation theory, namely,

$$\begin{aligned} \delta _{1} = 0, \quad \delta _{2} = 0, \end{aligned}$$
(33)

that makes \(R^{(\ell )}_{\mathrm{appr}}(s)\) (31) identical to the naive expression \(R^{(\ell )}_{\mathrm{naive}}(s)\) (12) for \(\ell =1\) and \(\ell =2\). Starting from the third order of perturbation theory (i.e., for \(\ell \ge 3\)) the coefficients \(\delta _j\) (32) are no longer vanishing and constitute a combination of the pertinent perturbative expansion coefficients \(d_{j}\) and \(\beta _{j}\) of the first \((\ell -2)\) orders. Specifically, the third-order and fourth-order coefficients read [42, 118, 119]

$$\begin{aligned} \delta _{3} = \frac{\pi ^2}{3} d_{1} B_{0}^{2} = \frac{\pi ^2}{3} d_{1} \end{aligned}$$
(34)

and

$$\begin{aligned} \delta _{4} = \frac{\pi ^2}{3} \left( \frac{5}{2} d_{1} B_{0}B_{1} + 3 d_{2}B_{0}^{2} \right) = \frac{\pi ^2}{3} \left( \frac{5}{2} d_{1} B_{1} + 3 d_{2} \right) ,\nonumber \\ \end{aligned}$$
(35)

respectively. In these equations \(d_j\) denote the Adler function perturbative expansion coefficients (8), whereas \(B_{j}=\beta _j/\beta ^{j+1}_{0}\) stands for the combination of perturbative coefficients of the renormalization group \(\beta \) function. In turn, at the fifth and sixth orders the coefficients \(\delta _j\) (32) can be represented as [31, 42, 121]

$$\begin{aligned} \delta _{5}= & {} \frac{\pi ^2}{3} \left[ \frac{3}{2} d_{1} ( B_{1}^{2} + 2 B_{2})+ 7 d_{2} B_{1} + 6 d_3\right] - \frac{\pi ^4}{5} d_{1},\nonumber \\ \end{aligned}$$
(36)
$$\begin{aligned} \delta _{6}= & {} \frac{\pi ^2}{3} \left[ \frac{7}{2} d_{1} (B_{1} B_{2} + B_{3})+4 d_{2}(B_{1}^{2} + 2 B_{2})\right. \nonumber \\&\left. +\,\frac{27}{2} d_3 B_{1} + 10 d_{4} \right] - \frac{\pi ^4}{5} \left( \frac{77}{12} d_{1} B_{1} + 5 d_{2}\right) . \end{aligned}$$
(37)
Table 2 Numerical values of the coefficients \(\delta _{j}\) (32) embodying the contributions of the relevant \(\pi ^2\)-terms (30). The last column employs the numerical estimation of the Adler function perturbative expansion coefficient \(d_5\) [42]
Table 3 Numerical values of the coefficients \(r_{j}\) of the re-expanded approximate R-ratio (31). The last column employs the numerical estimation of the Adler function perturbative expansion coefficient \(d_5\) [42]

The seventh- and eighth-order coefficients \(\delta _j\) (32) have recently been calculated as well, specifically [31, 121]

$$\begin{aligned} \delta _{7}= & {} \frac{\pi ^2}{3} \left[ 4 d_{1} \left( B_{1} B_{3} + \frac{1}{2} B_{2}^{2} + B_{4} \right) + 9 d_{2} (B_{1} B_{2} + B_{3})\right. \nonumber \\&\left. +\, \frac{15}{2} d_{3} (B_{1}^{2} + 2 B_{2})+ 22 d_{4} B_{1} + 15 d_{5} \right] \nonumber \\&-\, \frac{\pi ^4}{5} \left[ \frac{5}{6} d_{1} (17 B_{1}^{2} + 12 B_{2}) + \frac{57}{2} d_{2} B_{1} + 15 d_{3} \right] \nonumber \\&+\,\frac{\pi ^6}{7} d_1, \end{aligned}$$
(38)
$$\begin{aligned} \delta _{8}= & {} \frac{\pi ^2}{3} \left[ \frac{9}{2} d_{1} (B_{1} B_{4} + B_{2} B_{3} + B_{5})\right. \nonumber \\&+\, 10 d_{2} \left( B_{1} B_{3} + \frac{1}{2} B_{2}^{2} + B_{4} \right) + \frac{33}{2} d_{3}(B_{1} B_{2} + B_{3})\nonumber \\&\left. +\, 12 d_{4} (B_{1}^{2} + 2 B_{2})+\frac{65}{2} d_{5} B_{1} + 21 d_{6}\right] \nonumber \\&-\, \frac{\pi ^4}{5} \left[ \frac{15}{8} d_{1} (7 B_{1}^{3} + 22 B_{1} B_{2} + 8 B_{3})\right. \nonumber \\&\left. +\,\frac{5}{12} d_{2}(139 B_{1}^{2} + 96 B_{2})+ \frac{319}{4} d_{3} B_{1} + 35 d_{4}\right] \nonumber \\&+\,\frac{\pi ^6}{7} \left( \frac{223}{20} d_{1} B_{1} + 7 d_{2} \right) . \end{aligned}$$
(39)

The explicit expressions for the coefficients \(\delta _j\) (32) at the higher orders can be found in App. C of Ref. [31].

Table 4 The relative weight \((1+|d_{j}/\delta _{j}|)^{-1} \times 100\%\) of the \(\pi ^2\)-terms in the coefficients \(r_j\) of the re-expanded approximate R-ratio (31). The last column employs the numerical estimation of the Adler function perturbative expansion coefficient \(d_5\) [42]

Table 2 presents the numerical values of the first seven coefficients \(\delta _j\) (32), which embody the contributions of the relevant \(\pi ^2\)-terms (30). Tables 1 and 2 make it evident that in Eq. (31) the coefficients \(\delta _j\) can in no way be regarded as small corrections to the Adler function perturbative expansion coefficients \(d_j\) (8) for \(j \ge 3\). On the contrary, the values of coefficients \(\delta _j\) significantly exceed the values of respective perturbative coefficients \(d_j\), thereby constituting the dominant contribution to the coefficients \(r_j\) (32); see also Tables 3 and 4. As will be discussed below, eventually this results in an essential distortion of the re-expanded approximation \(R^{(\ell )}_{\text {appr}}(s)\) (31) with respect to the naive expression \(R^{(\ell )}_{\text {naive}}(s)\) (12). In particular, the higher-order terms of \(R^{(\ell )}_{\text {appr}}(s)\) (31) turn out to be substantially amplified and even sign-reversed with respect to those of \(R^{(\ell )}_{\text {naive}}(s)\) (12); see Tables 1 and 3. It is worthwhile to note also that, as one can infer from Tables 2 and 4, the values of coefficients \(\delta _j\) rapidly increase as the order j increases, which makes the loop convergence of \(R^{(\ell )}_{\text {appr}}(s)\) (31) worse than that of both \(R^{(\ell )}(s)\) (10) and \(R^{(\ell )}_{\text {naive}}(s)\) (12); see Sect. 6 for the details.

5 Spectral function at the higher-loop levels

It is certainly desirable to leave the truncated re-expanded approximation \(R^{(\ell )}_{\text {appr}}(s)\) (31) aside and have the function \(R^{(\ell )}(s)\) (10) calculated in a straightforward way beyond the one-loop level. To achieve this objective, the corresponding explicit expression for the involved spectral function \(\rho ^{(\ell )}(\sigma )\) (11) is required. Despite the latter becomes rather cumbrous for \(\ell \ge 2\), the following method enables one to calculate \(\rho ^{(\ell )}(\sigma )\) explicitly at an arbitraryFootnote 6 loop level.

Specifically, it proves to be convenient to express the spectral function \(\rho ^{(\ell )}(\sigma )\) (11) in terms of the so-called “partial” spectral functions \({\bar{\rho }}^{(\ell )}_{j}(\sigma )\) corresponding to the jth power of the \(\ell \)-loop perturbative couplant (9), namely

$$\begin{aligned} \rho ^{(\ell )}(\sigma ) = \sum _{j=1}^{\ell } d_{j}\, {\bar{\rho }}^{(\ell )}_{j}(\sigma ), \end{aligned}$$
(40)

where

$$\begin{aligned}&{\bar{\rho }}^{(\ell )}_{j}(\sigma ) \nonumber \\&\quad = \frac{1}{2\pi i} \lim _{\varepsilon \rightarrow 0_{+}} \left\{ \left[ a^{(\ell )}_{\text {s}} ({-}\sigma {-}i\varepsilon )\right] ^{j} {-} \left[ a^{(\ell )}_{\text {s}} ({-}\sigma {+}i\varepsilon )\right] ^{j} \right\} .\nonumber \\ \end{aligned}$$
(41)

To obtain the partial spectral function \({\bar{\rho }}^{(\ell )}_{j}(\sigma )\) (41) for any \(j \ge 1\) it appears to be enough to calculate only the real and imaginary parts of the \(\ell \)-loop couplant \(a^{(\ell )}_{\text {s}}(Q^2)\) at the edges of its cut:

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0_{+}} a^{(\ell )}_{\text {s}}(-\sigma \pm i\varepsilon ) = a^{(\ell )}_{\text { Re}}(\sigma ) \mp i \pi a^{(\ell )}_{\text { Im}}(\sigma ). \end{aligned}$$
(42)

Here \(a^{(\ell )}_{\text { Re}}(\sigma )\) and \(a^{(\ell )}_{\text { Im}}(\sigma )\) are the real functions of their arguments and \(\sigma \ge 0\) is assumed. For positive integer values of j the following equation holds:

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0_{+}} \Bigl [a^{(\ell )}_{\text {s}}(-\sigma \pm i\varepsilon )\Bigr ]^{j}= & {} \sum _{k=0}^{j} \left( {\begin{array}{c}j\\ k\end{array}}\right) \, (\mp i\pi )^{k}\nonumber \\&\times \Bigl [a^{(\ell )}_{\text { Re}}(\sigma )\Bigr ]^{j-k}\, \Bigl [a^{(\ell )}_{\text { Im}}(\sigma )\Bigr ]^{k}, \end{aligned}$$
(43)

with

$$\begin{aligned} \left( {\begin{array}{c}n\\ m\end{array}}\right) = \frac{n!}{m!\,(n-m)!} \end{aligned}$$
(44)

being the binomial coefficient. Then it is convenient to isolate in Eq. (43) its real and imaginary parts, specifically

$$\begin{aligned}&\lim _{\varepsilon \rightarrow 0_{+}} \Bigl [a^{(\ell )}_{\text {s}}(-\sigma \pm i\varepsilon )\Bigr ]^{j}=\sum _{k=0}^{K(j+1)} \left( {\begin{array}{c}j\\ 2k\end{array}}\right) \, (-1)^{k}(\pi )^{2k}\nonumber \\&\quad \times \, \Bigl [a^{(\ell )}_{\text { Re}}(\sigma )\Bigr ]^{j-2k}\, \Bigl [a^{(\ell )}_{\text { Im}}(\sigma )\Bigr ]^{2k} \nonumber \\&\quad \mp \, i \pi \sum _{k=0}^{K(j)} \left( {\begin{array}{c}j\\ 2k+1\end{array}}\right) \, (-1)^{k}(\pi )^{2k} \nonumber \\&\quad \times \, \Bigl [a^{(\ell )}_{\text { Re}}(\sigma )\Bigr ]^{j-2k-1}\, \Bigl [a^{(\ell )}_{\text { Im}}(\sigma )\Bigr ]^{2k+1}, \end{aligned}$$
(45)

where

$$\begin{aligned} K(j) = \frac{j-2}{2} + \frac{j \;\hbox {mod}\; 2}{2} \end{aligned}$$
(46)

and \((j \;\hbox {mod}\; n)\) denotes the remainder on division of j by n. Therefore, the partial spectral function (41) reads (see also Refs. [107, 122, 123])

$$\begin{aligned} {\bar{\rho }}^{(\ell )}_{j}(\sigma )= & {} \sum _{k=0}^{K(j)} \left( {\begin{array}{c}j\\ 2k+1\end{array}}\right) \, (-1)^{k}\, \pi ^{2k}\nonumber \\&\times \,\Bigl [a^{(\ell )}_{\text { Re}}(\sigma )\Bigr ]^{j-2k-1}\, \Bigl [a^{(\ell )}_{\text { Im}}(\sigma )\Bigr ]^{2k+1}, \end{aligned}$$
(47)

with \(\sigma \ge 0\) and \(j \ge 1\) being assumed. In particular, the first five relations (47) acquire a compact form

$$\begin{aligned} {\bar{\rho }}^{(\ell )}_{1}(\sigma )= & {} a^{(\ell )}_{\text { Im}}(\sigma ), \end{aligned}$$
(48)
$$\begin{aligned} {\bar{\rho }}^{(\ell )}_{2}(\sigma )= & {} 2\, a^{(\ell )}_{\text { Im}}(\sigma )\, a^{(\ell )}_{\text { Re}}(\sigma ), \end{aligned}$$
(49)
$$\begin{aligned} {\bar{\rho }}^{(\ell )}_{3}(\sigma )= & {} a^{(\ell )}_{\text { Im}}(\sigma ) \left\{ 3\Bigl [a^{(\ell )}_{\text { Re}}(\sigma )\Bigr ]^{2} - \pi ^{2}\Bigl [a^{(\ell )}_{\text { Im}}(\sigma )\Bigr ]^{2}\right\} , \end{aligned}$$
(50)
$$\begin{aligned} {\bar{\rho }}^{(\ell )}_{4}(\sigma )= & {} 4\, a^{(\ell )}_{\text { Im}}(\sigma )\, a^{(\ell )}_{\text { Re}}(\sigma ) \left\{ \Bigl [a^{(\ell )}_{\text { Re}}(\sigma )\Bigr ]^{2} {-} \pi ^{2}\Bigl [a^{(\ell )}_{\text { Im}}(\sigma )\Bigr ]^{2}\right\} ,\nonumber \\ \end{aligned}$$
(51)
$$\begin{aligned} {\bar{\rho }}^{(\ell )}_{5}(\sigma )= & {} a^{(\ell )}_{\text { Im}}(\sigma ) \biggl \{ 5\Bigl [a^{(\ell )}_{\text { Re}}(\sigma )\Bigr ]^{4} - 10 \pi ^{2}\Bigl [a^{(\ell )}_{\text { Im}}(\sigma )a^{(\ell )}_{\text { Re}}(\sigma )\Bigr ]^{2} \nonumber \\&+ \pi ^{4}\Bigl [a^{(\ell )}_{\text { Im}}(\sigma )\Bigr ]^{4}\biggr \}. \end{aligned}$$
(52)

In turn, the \(\ell \)-loop functions \(a^{(\ell )}_{\text { Re}}(\sigma )\) and \(a^{(\ell )}_{\text { Im}}(\sigma )\) entering Eq. (47) can also be explicitly calculated in a similar way. Specifically, as mentioned in Sect. 2, the perturbative QCD couplant \(a^{(\ell )}_{\text {s}}(Q^2)\) can be represented as the double sum (9) comprised of the functions

$$\begin{aligned} {\bar{a}}_{n}^{m}(Q^2) = \frac{\ln ^{m}(\ln z)}{\ln ^n z}, \quad z = \frac{Q^2}{\varLambda ^2}. \end{aligned}$$
(53)

It is worthwhile to decompose the function \(\bar{a}_{n}^{m}(Q^2)\) (53) at the edges of its cut into the real and imaginary parts, namely

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0_{+}} {\bar{a}}_{n}^{m}(-\sigma \pm i\varepsilon ) = u_{n}^{m}(\sigma ) \mp i \pi v_{n}^{m}(\sigma ), \end{aligned}$$
(54)

where \(u_{n}^{m}(\sigma )\) and \(v_{n}^{m}(\sigma )\) are the real functions of their arguments and \(\sigma \ge 0\) is assumed. Therefore, the functions \(a^{(\ell )}_{\text { Re}}(\sigma )\) and \(a^{(\ell )}_{\text { Im}}(\sigma )\) (42) take the following form:

$$\begin{aligned} a^{(\ell )}_{\text { Re}}(\sigma )= & {} \sum _{n=1}^{\ell }\sum _{m=0}^{n-1} b^{m}_{n}\, u_{n}^{m}(\sigma ), \end{aligned}$$
(55)
$$\begin{aligned} a^{(\ell )}_{\text { Im}}(\sigma )= & {} \sum _{n=1}^{\ell }\sum _{m=0}^{n-1} b^{m}_{n}\, v_{n}^{m}(\sigma ). \end{aligned}$$
(56)

On the left-hand side of Eq. (53) and in Eqs. (54)–(56) the integer superscripts m are not to be confused with respective powers, whereas the coefficients \(b^{m}_{n}\) have been specified in Eq. (9).

Then, to calculate the functions \(u_{n}^{m}(\sigma )\) and \(v_{n}^{m}(\sigma )\) entering Eqs. (55) and (56), it is convenient to split the left-hand side of Eq. (54) into two factors:

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0_{+}} {\bar{a}}_{n}^{m}(-\sigma \pm i\varepsilon ) = \lim _{\varepsilon \rightarrow 0_{+}} \Bigl [{\bar{a}}_{n}^{0}(-\sigma \pm i\varepsilon )\, {\bar{a}}_{0}^{m}(-\sigma \pm i\varepsilon ) \Bigr ].\nonumber \\ \end{aligned}$$
(57)

The first factor on the right-hand side of Eq. (57) reads

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0_{+}} {\bar{a}}_{n}^{0}(-\sigma \pm i\varepsilon ) = \frac{(y \mp i\pi )^{n}}{(y^2 + \pi ^2)^{n}}, \quad y=\ln \left( \frac{\sigma }{\varLambda ^2}\right) . \end{aligned}$$
(58)

Proceeding along the same lines as earlier, one can cast the numerator on the right-hand side of this equation into

$$\begin{aligned} (y \mp i\pi )^{n}= & {} \sum _{k=0}^{K(n+1)} \left( {\begin{array}{c}n\\ 2k\end{array}}\right) (-1)^{k} \pi ^{2k} y^{n-2k} \nonumber \\&\mp \,i\pi \sum _{k=0}^{K(n)} \left( {\begin{array}{c}n\\ 2k+1\end{array}}\right) (-1)^{k} \pi ^{2k} y^{n-2k-1}, \end{aligned}$$
(59)

with K(n) being specified in Eq. (46). Hence, the functions \(u_{n}^{m}(\sigma )\) and \(v_{n}^{m}(\sigma )\) (54) for \(m=0\) read

$$\begin{aligned} u_{n}^{0}(\sigma )= & {} \frac{1}{(y^2 + \pi ^2)^{n}} \sum _{k=0}^{K(n+1)} \left( {\begin{array}{c}n\\ 2k\end{array}}\right) (-1)^{k} \pi ^{2k} y^{n-2k}, \end{aligned}$$
(60)
$$\begin{aligned} v_{n}^{0}(\sigma )= & {} \frac{1}{(y^2 + \pi ^2)^{n}} \sum _{k=0}^{K(n)} \left( {\begin{array}{c}n\\ 2k+1\end{array}}\right) (-1)^{k} \pi ^{2k} y^{n-2k-1},\nonumber \\ \end{aligned}$$
(61)

with \(n \ge 1\) being assumed.

In turn, for \(m \ge 1\) the second factor on the right-hand side of Eq. (57) takes the following form:

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0_{+}} {\bar{a}}_{0}^{m}(-\sigma \pm i\varepsilon ) = \Bigl [\ln (y \pm i\pi )\Bigr ]^{m}, \end{aligned}$$
(62)

where \(y=\ln (\sigma /\varLambda ^2)\). Since for real a and b

$$\begin{aligned}&\ln (a \pm i b) = \ln \sqrt{ a^2 + b^2} \pm i\pi \left[ \frac{1}{2} - \frac{1}{\pi }\arctan \left( \frac{a}{b}\right) \right] ,\nonumber \\&\quad b>0, \end{aligned}$$
(63)

Eq. (62) can be rewritten as

$$\begin{aligned}&\lim _{\varepsilon \rightarrow 0_{+}} {\bar{a}}_{0}^{m}(-\sigma \pm i\varepsilon ) \nonumber \\&\quad = \sum \limits _{k=0}^{m}\left( {\begin{array}{c}m\\ k\end{array}}\right) (\pm i\pi )^{k} [L_{1}(y)]^{m-k}\, [L_{2}(y)]^{k}, \end{aligned}$$
(64)

where

$$\begin{aligned} L_{1}(y) = \ln \sqrt{y^{2}+\pi ^{2}}, \quad L_{2}(y) = \frac{1}{2} - \frac{1}{\pi }\arctan \left( \frac{y}{\pi }\right) . \end{aligned}$$
(65)

Following the same steps as above, one can cast Eq. (64) into

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0_{+}} {\bar{a}}_{0}^{m}(-\sigma \pm i\varepsilon ) = u_{0}^{m}(\sigma ) \mp i \pi v_{0}^{m}(\sigma ), \end{aligned}$$
(66)

where

$$\begin{aligned} u_{0}^{m}(\sigma )= & {} \sum \limits _{k=0}^{K(m+1)}\left( {\begin{array}{c}m\\ 2k\end{array}}\right) (-1)^{k}\pi ^{2k} \nonumber \\&\times \, [L_{1}(y)]^{m-2k}\, [L_{2}(y)]^{2k}, \end{aligned}$$
(67)
$$\begin{aligned} v_{0}^{m}(\sigma )= & {} \sum \limits _{k=0}^{K(m)}\left( {\begin{array}{c}m\\ 2k+1\end{array}}\right) (-1)^{k+1}\pi ^{2k}\nonumber \\&\times \, [L_{1}(y)]^{m-2k-1}\,[L_{2}(y)]^{2k+1}, \end{aligned}$$
(68)

K(m) is defined in Eq. (46), and \(m \ge 1\) is assumed.

Therefore, the functions \(u_{n}^{m}(\sigma )\) and \(v_{n}^{m}(\sigma )\) (54) read

$$\begin{aligned} u_{n}^{m}(\sigma ) = {\left\{ \begin{array}{ll} u_{n}^{0}(\sigma ),&{} \hbox {if } m=0,\\ u_{n}^{0}(\sigma ) u_{0}^{m}(\sigma ) - \pi ^2 v_{n}^{0}(\sigma ) v_{0}^{m}(\sigma ),&{}\hbox {if }m \ge 1, \end{array}\right. } \end{aligned}$$
(69)

and

$$\begin{aligned} v_{n}^{m}(\sigma ) = {\left\{ \begin{array}{ll} v_{n}^{0}(\sigma ),&{} \hbox {if }m=0,\\ v_{n}^{0}(\sigma ) u_{0}^{m}(\sigma ) + u_{n}^{0}(\sigma ) v_{0}^{m}(\sigma ),&{} \hbox {if }m \ge 1, \end{array}\right. } \end{aligned}$$
(70)

with \(n \ge 1\) being assumed. Thus, the explicit expression for the \(\ell \)-loop spectral function \(\rho ^{(\ell )}(\sigma )\) (40) takes the following form:

$$\begin{aligned} \rho ^{(\ell )}(\sigma )= & {} \sum _{j=1}^{\ell } d_{j} \sum _{k=0}^{K(j)} \left( {\begin{array}{c}j\\ 2k+1\end{array}}\right) (-1)^{k}\, \pi ^{2k}\nonumber \\&\times \left[ \sum _{n=1}^{\ell }\sum _{m=0}^{n-1} b^{m}_{n}\, u_{n}^{m}(\sigma )\right] ^{j-2k-1}\nonumber \\&\times \left[ \sum _{n=1}^{\ell }\sum _{m=0}^{n-1} b^{m}_{n}\, v_{n}^{m}(\sigma )\right] ^{2k+1}, \end{aligned}$$
(71)

where the functions \(u_{n}^{m}(\sigma )\) and \(v_{n}^{m}(\sigma )\) are specified in Eqs. (69) and (70), respectively.

Fig. 4
figure 4

The partial spectral functions \({\bar{\rho }}^{(\ell )}_{j}(\sigma )\) (47) corresponding to the j-th power (\(1 \le j \le \ell \)) of the \(\ell \)-loop (\(2 \le \ell \le 5\)) perturbative QCD couplant \(a^{(\ell )}_{\mathrm{s}}(Q^2)\). Plot A: two-loop level (\(\ell =2\), \(1 \le j \le 2\)). Plot B: three-loop level (\(\ell =3\), \(1 \le j \le 3\)), the function \({\bar{\rho }}^{(3)}_{3}(\sigma )\) is scaled by the factor of 10. Plot C: four-loop level (\(\ell =4\), \(1 \le j \le 4\)), the functions \({\bar{\rho }}^{(4)}_{3}(\sigma )\) and \({\bar{\rho }}^{(4)}_{4}(\sigma )\) are scaled by the factors of 10 and \(10^2\), respectively. Plot D: five-loop level (\(\ell =5\), \(1 \le j \le 5\)), the functions \({\bar{\rho }}^{(5)}_{3}(\sigma )\), \({\bar{\rho }}^{(5)}_{4}(\sigma )\), and \({\bar{\rho }}^{(5)}_{5}(\sigma )\) are scaled by the factors of 10, \(10^2\), and \(10^2\), respectively

The higher-loop partial spectral functions \({\bar{\rho }}^{(\ell )}_{j}(\sigma )\) (47), which correspond to the jth power (\(1 \le j \le \ell \)) of the \(\ell \)-loop (\(2 \le \ell \le 5\)) perturbative QCD couplant \(a^{(\ell )}_{\text {s}}(Q^2)\), are displayed in Fig. 4. As one can infer from this figure, the functions \({\bar{\rho }}^{(\ell )}_{j}(\sigma )\) (47) vanish at both \(\sigma \rightarrow \infty \) and \(\sigma \rightarrow 0\), the higher-order functions \({\bar{\rho }}^{(\ell )}_{j}(\sigma )\) (\(j \ge 2\)) being substantially suppressed with respect to those of the preceding orders. For example, Fig. 4D implies that at the five-loop level the maximum value of the fifth-order function \({\bar{\rho }}^{(5)}_{5}(\sigma )\) is about three orders of magnitude less than the maximum value of the first-order function \({\bar{\rho }}^{(5)}_{1}(\sigma )\). In turn, as it will be discussed in the next section, the fact that the function \({\bar{\rho }}^{(\ell )}_{j+1}(\sigma )\) is subdominant to \({\bar{\rho }}^{(\ell )}_{j}(\sigma )\) eventually results in an enhanced higher-loop stability of the proper expression for the R-ratio (10) at moderate and low energies with respect to both its naive form (12) and the commonly employed truncated re-expanded approximation (31).

Fig. 5
figure 5

Plot A: The spectral function \(\rho ^{(\ell )}(\sigma )\) (71) at the first five loop levels (\(1 \le \ell \le 5\)). Plot B: The relative difference \(\rho ^{(\ell )}_{\text {diff}}(\sigma )\) (74) between the \(\ell \)-loop and \((\ell +1)\)-loop spectral functions (71) at various loop levels. The function \(\rho ^{(4)}_{\text {diff}}(\sigma )\) is scaled by the factor of 10

The plots of the spectral function \(\rho ^{(\ell )}(\sigma )\) (71) at the first five loop levels (\(1 \le \ell \le 5\)) are displayed in Fig. 5A. As one can infer from this figure, the function \(\rho ^{(\ell )}(\sigma )\) vanishes at both \(\sigma \rightarrow \infty \) and \(\sigma \rightarrow 0\). Specifically, at the higher-loop levels (\(\ell \ge 2\))

$$\begin{aligned} \rho ^{(\ell )}(\sigma ) \simeq \frac{d_{1}}{y^2} + \mathcal {O}\left( \frac{1}{y^3}\right) , \quad y=\ln \left( \frac{\sigma }{\varLambda ^2}\right) , \quad y\rightarrow \infty \nonumber \\ \end{aligned}$$
(72)

and

$$\begin{aligned} \rho ^{(\ell )}(\sigma ) \simeq \frac{d_{1}(1+B_{1})}{y^2} + \mathcal {O}\left( \frac{1}{y^3}\right) , \qquad y \rightarrow -\infty . \end{aligned}$$
(73)

Additionally, Fig. 5A implies that the range of y, where the difference between \(\rho ^{(\ell )}(\sigma )\) and \(\rho ^{(\ell +1)}(\sigma )\) is sizable, is located in the vicinity of \(y=0\) and becomes smaller at larger \(\ell \). This issue is also elucidated by Fig. 5B, which shows the relative difference between the \(\ell \)-loop and \((\ell +1)\)-loop spectral functions (71)

$$\begin{aligned} \rho ^{(\ell )}_{\text {diff}}(\sigma ) = \left[ 1-\frac{\rho ^{(\ell )}(\sigma )}{\rho ^{(\ell +1)}(\sigma )}\right] \times 100\% \end{aligned}$$
(74)

at various loop levels.

6 R-ratio at the higher-loop levels

The obtained in the previous section explicit expression for the \(\ell \)-loop spectral function \(\rho ^{(\ell )}(\sigma )\) (71) enables one to calculate the R-ratio of electron–positron annihilation into hadrons (10) at an arbitrary loop level (assuming that the involved perturbative coefficients \(\beta _j\) and \(d_j\) are available). The integration in Eq. (10) can be directly performed for the function \(\rho ^{(\ell )}(\sigma )\) (71) as a whole, though for the illustrative purposes it is somewhat convenient to keep the contributions of the partial spectral functions \({\bar{\rho }}^{(\ell )}_{j}(\sigma )\) (47) to Eq. (40) separately from each other, which casts Eq. (10) into

$$\begin{aligned} R^{(\ell )}(s) = 1 + r^{(\ell )}(s), \quad r^{(\ell )}(s) = \sum \limits _{j=1}^{\ell } d_{j}\,A^{(\ell )}_{\text { TL},j}(s). \end{aligned}$$
(75)

In this equation

$$\begin{aligned} A^{(\ell )}_{\text { TL},j}(s) = \int \limits _{s}^{\infty } {\bar{\rho }}^{(\ell )}_{j}(\sigma )\, \frac{\mathrm{d}\sigma }{\sigma } \end{aligned}$$
(76)

stands for the jth-order \(\ell \)-loop “timelike” effective couplant, which constitutes the proper continuation of the jth power of \(\ell \)-loop QCD couplant \([a^{(\ell )}_{\text {s}}(Q^2)]^{j}\) into the timelike domain.

As discussed in Sect. 3, at the one-loop level the first-order function (76) acquires a quite simple form, namely (17)

$$\begin{aligned} A^{(1)}_{\text { TL},1}(s) = a^{(1)}_{\text { TL}}(s) = \frac{1}{2} - \frac{1}{\pi }\arctan \left( \frac{\ln w}{\pi }\right) , \quad w = \frac{s}{\varLambda ^2}.\nonumber \\ \end{aligned}$$
(77)

Despite the fact that the partial spectral functions \({\bar{\rho }}^{(\ell )}_{j}(\sigma )\) (47) become rather cumbrous at the higher loop levels, it appears that the integration in Eq. (76) can be performed explicitly for \(\ell >1\), too. For example, the two-loop first-order function (76) reads [9, 10]

$$\begin{aligned}&A^{(2)}_{\text { TL},1}(s) = a^{(1)}_{\text { TL}}(s) - \frac{B_{1}}{\ln ^{2}w+\pi ^2}\nonumber \\&\quad \times \left[ W(s) - a^{(1)}_{\text { TL}}(s) \ln w + 1 \right] ,\nonumber \\ \end{aligned}$$
(78)

whereas the second-order one can be represented as

$$\begin{aligned} A^{(2)}_{\text { TL},2}(s)= & {} \frac{1}{\ln ^{2}w+\pi ^2} + \frac{B_{1}}{\left( \ln ^{2}w+\pi ^2\right) ^2} \nonumber \\&\times \left\{ a^{(1)}_{\text { TL}}(s)\Bigl (\ln ^{2}w-\pi ^2\Bigr ) - \ln w\,[2W(s)+1]\right\} \nonumber \\&+\,\frac{B_{1}^{2}}{(\ln ^{2}w+\pi ^2)^{ 3}} \left\{ \left( \ln ^{2}w - \frac{\pi ^2}{3}\right) \right. \nonumber \\&\times \left[ \left( W(s)+\frac{1}{3}\right) ^{ 2} + \frac{1}{9} - \pi ^2\left[ a^{(1)}_{\text { TL}}(s)\right] ^{2}\right] \nonumber \\&\left. -\,\frac{2}{3}a^{(1)}_{\text { TL}}(s)\ln w(\ln ^{2}w {-} 3\pi ^2) \left( W(s){+}\frac{1}{3}\right) \right\} .\nonumber \\ \end{aligned}$$
(79)

In these equations \(a^{(1)}_{\text { TL}}(s)\) is given by Eq. (77), \(B_j=\beta _j/\beta ^{j+1}_{0}\), and

$$\begin{aligned} W(s) = \ln \sqrt{\ln ^{2}w + \pi ^2},\quad w = \frac{s}{\varLambda ^2}. \end{aligned}$$
(80)

In turn, at the three-loop level the first-order function (76) takes the form

$$\begin{aligned} A^{(3)}_{\text { TL},1}(s)= & {} A^{(2)}_{\text { TL},1}(s) - \frac{B_{1}^2}{\left( \ln ^{2}w+\pi ^2\right) ^{ 2}} [T_{1}(s)T_{2}(s) + \ln w]\nonumber \\&+\frac{B_{2}\ln w}{(\ln ^{2}w+\pi ^2)^{ 2}}, \end{aligned}$$
(81)

where \(A^{(2)}_{\text { TL},1}(s)\) is specified in Eq. (78) and

$$\begin{aligned} T_{1}(s)= & {} a^{(1)}_{\text { TL}}(s)\,\ln w - W(s), \end{aligned}$$
(82)
$$\begin{aligned} T_{2}(s)= & {} \pi ^2 a^{(1)}_{\text { TL}}(s) + W(s)\,\ln w. \end{aligned}$$
(83)

As for the four-loop first-order function (76), it can be represented as

$$\begin{aligned} A^{(4)}_{\text { TL},1}(s)= & {} A^{(3)}_{\text { TL},1}(s) + \frac{1}{\left( \ln ^{2}w+\pi ^2\right) ^{ 3}}\frac{B_{3}}{2} \left( \ln ^{2}w-\frac{\pi ^2}{3}\right) \nonumber \\&+\,\frac{B_{1}B_{2}}{(\ln ^{2}w+\pi ^2)^{ 3}} \biggl \{\frac{\pi ^2}{3} - 3T_{2}(s)\ln w +\pi ^2 W(s) \nonumber \\&+\, \ln ^{2}w\Bigl [a^{(1)}_{\text { TL}}(s)\ln w-1 \Bigr ] \biggr \} - \frac{1}{\left( \ln ^{2}w+\pi ^2\right) ^{ 3}}\frac{B_{1}^3}{2} \nonumber \\&\times \Biggl \{ 2T_{2}(s)\ln w \Bigl [W(s)-3\Bigr ] -T_{2}^2(s) \Bigl [2T_{1}(s)+1\Bigr ]\nonumber \\&+\,\pi ^2T_{1}^2(s)\biggl [\frac{2}{3}T_{1}(s)+3\biggr ] + \ln ^{2}w\Bigl [2\,a^{(1)}_{\text { TL}}(s)\ln w-1\Bigr ]\nonumber \\&+2W^{2}(s)\ln ^{2}w\Bigl [T_{1}(s)+W(s)-2\Bigr ] + \frac{\pi ^2}{3}\nonumber \\&+2W(s)\Bigl [1-W(s)\Bigr ] \Bigl [\pi ^2+a^{(1)}_{\text { TL}}(s)\ln ^{3}w\Bigr ] \Biggr \}, \end{aligned}$$
(84)

where the functions \(A^{(3)}_{\text { TL},1}(s)\), W(s), \(T_{1}(s)\), and \(T_{2}(s)\) are given by Eqs. (81), (80), (82), and (83), respectively. It is worthwhile to note also that the functions \(A^{(\ell )}_{\text { TL},j}(s)\) (76) entering Eq. (75) can be computed numerically by making use of the routines included in the freely available program packages [122, 123] [which, being based on a less universal method of calculation of the relevant spectral function (11) than that of Eq. (71), is applicable at first four loop levels only] and [124].

Fig. 6
figure 6

Two-loop timelike effective expansion function of the second order \(A^{(2)}_{\text { TL},2}(s)\) [Eq. (79), solid curve] and the approximations corresponding to various orders of its re-expansion (85). The result of “naive” continuation (12) [the first term on the right-hand side of Eq. (85)] is shown by dashed curve. Numerical label next to a dot-dashed curve indicates the highest absolute value of power of \(\ln w\) retained on the right-hand side of Eq. (85). The plotted functions are scaled by the factor of 10. Vertical solid line marks the boundary of convergence range of the re-expansion (85) at \(\sqrt{s}/\varLambda = \exp (\pi /2) \simeq 4.81\)

As discussed earlier, the \(\pi ^2\)-terms play a valuable role in the studies of the strong interaction processes in the timelike domain, and their ignorance (complete or partial) may yield misleading results. In particular, this issue is illustrated by Fig. 6, which displays the two-loop timelike effective expansion function of the second order \(A^{(2)}_{\text { TL},2}(s)\) [Eq. (79), solid curve] and the approximations (dashed and dot-dashed curves) corresponding to various orders of its re-expansion for \(\sqrt{s}/\varLambda > \exp (\pi /2) \simeq 4.81\):

$$\begin{aligned} A^{(2)}_{\text { TL},2}(s)\simeq & {} \left[ a^{(2)}_{\text {s}}(|s|)\right] ^{2} - \frac{\pi ^2}{\ln ^{4}w} + \frac{\pi ^2}{\ln ^{5}w}\,B_{1} \left( 4\ln \ln w-\frac{7}{3}\right) \nonumber \\&-\, \frac{\pi ^2}{\ln ^{6}w} \left[ \frac{B_{1}^{2}}{3}( 10\ln ^{2}\ln w - 9\ln \ln w + 1 ) - \pi ^2\right] \nonumber \\&-\, \frac{\pi ^4}{\ln ^{7}w}\frac{3B_{1}}{10}(20\ln \ln w - 19) + \mathcal {O}\left( \frac{1}{\ln ^{8}w}\right) .\nonumber \\ \end{aligned}$$
(85)

Specifically, the dashed curve shows the result of naive continuation of the respective term of the Adler function perturbative expansion into the timelike domain (12) [the first term on the right-hand side of Eq. (85)], whereas the dot-dashed curves additionally includeFootnote 7 the lowest-order \(\pi ^2\)-terms [numerical labels indicate the highest absolute value of power of \(\ln w\) retained in Eq. (85)]. Figure 6 implies that the re-expansion (85) converges rather slowly at low and moderate energies. Furthermore, even at relatively high energies \(A^{(2)}_{\text { TL},2}(s)\) (79) considerably differs from \([a^{(2)}_{\text {s}}(|s|)]^{2}\), the latter being the only part of the function (79), which is retained in the two-loop approximation of the R-ratio (31). For example, as one can infer from Fig. 6, for \(\sqrt{s}/\varLambda =20\) the function \([a^{(2)}_{\text {s}}(|s|)]^{2}\) exceeds \(A^{(2)}_{\text { TL},2}(s)\) by about \(21\%\), and to securely achieve \(10\%\) accuracy in the re-expansion (85) the inclusion of the \(\pi ^2\)-terms up to the order of \(\ln ^{-7}w\) is required.

Fig. 7
figure 7

The function \(R^{(\ell )}(s)\) [Eq. (10), solid curves], its naive form [Eq. (12), dot-dashed curves], and the re-expanded approximation [Eq. (31), dashed curves] at various loop levels (\(1 \le \ell \le 5\)). Vertical solid line marks the boundary of convergence range of the re-expansion (31) at \(\sqrt{s}/\varLambda = \exp (\pi /2) \simeq 4.81\). Numerical labels specify the loop level

The issue of the higher-loop stability of the R-ratio of electron–positron annihilation into hadrons is also elucidated by Fig. 7. In particular, this figure displays the proper expression \(R^{(\ell )}(s)\) [Eq. (10), solid curves], its naive form \(R^{(\ell )}_{\text {naive}}(s)\) [Eq. (12), dot-dashed curves], and the truncated re-expanded approximation \(R^{(\ell )}_{\text {appr}}(s)\) [Eq. (31), dashed curves] at various loop levels (\(1 \le \ell \le 5\)). Figure 7 makes it evident that the loop convergence of the commonly employed approximation \(R^{(\ell )}_{\text {appr}}(s)\) (31) is worse than that of both the proper expression \(R^{(\ell )}(s)\) (10) and the naive one \(R^{(\ell )}_{\text {naive}}(s)\) (12). As discussed earlier, this is primarily caused by the fact that the convergence range of the re-expanded approximation (31) is strictly limited to \(\sqrt{s}/\varLambda > \exp (\pi /2) \simeq 4.81\), which, in turn, results in rather large values of the higher-order coefficients \(\delta _j\) (32) embodying the contributions of the corresponding \(\pi ^2\)-terms (30). In particular, as one can infer from Fig. 7, beyond the two-loop level (i.e., for \(\ell \ge 3\)) the curves corresponding to \(R^{(\ell )}(s)\) (10) are nearly indistinguishable from each other, whereas the curves corresponding to \(R^{(\ell )}_{\text {appr}}(s)\) (31) start to swerve quite above the boundary of its convergence range. For example, at \(\sqrt{s}/\varLambda = 2\exp (\pi /2) \simeq 9.62\) the relative difference between the \(\ell \)-loop and \((\ell +1)\)-loop strong corrections to the proper expression \(R^{(\ell )}(s)\) (10) is \(0.4\%\) for \(\ell =3\) and \(0.003\%\) for \(\ell =4\), whereas for the case of the commonly employed approximation \(R^{(\ell )}_{\text {appr}}(s)\) (31) these values increase up to \(4.4\%\) for \(\ell =3\) and \(2.5\%\) for \(\ell =4\).

Fig. 8
figure 8

The function \(R^{(\ell )}_{\text { diff}}(s)\) (86) in the energy range planned for the future CLIC experiment [4] at various loop levels

To expound the accuracy of approximation of the R-ratio (10) by its truncated re-expansion \(R^{(\ell )}_{\text {appr}}(s)\) (31), it is worthwhile to mention also the following. At a moderate energy scale of the mass of \(\tau \) lepton \(\sqrt{s}=M_{\tau }\) the relative difference between the \(\ell \)-loop strong corrections \(r^{(\ell )}(s)\) (10) and \(r^{(\ell )}_{\text {appr}}(s)\) (31) turns out to be as high as 26%, 28%, 14%, 2%, and 7% at the one-, two-, three-, four-, and five-loop levels, respectively. Moreover, it appears that even at high energies the ignorance of the higher-order \(\pi ^2\)-terms in the truncated re-expanded approximation (31) may produce a considerable effect. In particular, this issue is illustrated by Fig. 8, which displays the quantity

$$\begin{aligned} R^{(\ell )}_{\text {diff}}(s) = \left| \frac{R^{(\ell )}_{\text {appr}}(s)-R^{(\ell )}(s)}{R^{(\ell )}_{\text {appr}}(s)-R^{(\ell +1)}_{\text {appr}} (s)}\right| \times 100\% \end{aligned}$$
(86)

at various loop levels. Specifically, as one can infer from Fig. 8, in the energy range planned for the future CLIC experiment [4] the effect of inclusion of the \(\pi ^2\)-terms discarded in the approximate expression \(R^{(\ell )}_{\text {appr}}(s)\) (31) is either comparable to or prevailing over the effect of inclusion of the next-order perturbative correction.

7 Conclusions

The strong corrections to the R-ratio of electron–positron annihilation into hadrons are studied at the higher-loop levels. In particular, the derivation of a general form of the commonly employed approximate expression for the R-ratio (which constitutes its truncated re-expansion at high energies) is delineated, the appearance of the pertinent \(\pi ^2\)-terms is expounded, and their basic features are examined. It is demonstrated that the validity range of such approximation is strictly limited to \(\sqrt{s}/\varLambda > \exp (\pi /2) \simeq 4.81\) and that it converges rather slowly when the energy scale approaches this value. The spectral function required for the proper calculation of the R-ratio is explicitly derived and its properties at the higher-loop levels are studied. The developed method of calculation of the spectral function enables one to obtain the explicit expression for the latter at an arbitrary loop level. By making use of the derived spectral function the proper expression for the R-ratio is calculated up to the five-loop level and its properties are examined. In particular, it is shown that the loop convergence of the proper expression for the R-ratio is better than that of its commonly employed approximation and that the omitted higher-order \(\pi ^2\)-terms in the latter may produce a considerable effect.