1 Introduction

Quantum field theory (QFT) and the principle of gauge invariance stand as fundamental pillars of the Standard Model of particle physics, successfully describing three out of four forces in nature. QFT explains the forces in the macroscopic world as emerging from exchanges at a microscopic level of so-called virtual particles. These virtual particles mediate the interaction between real particles existing at macroscopic timescales. This concept of QFT allows an extremely precise theoretical description of phenomena seen in nature, e.g., running coupling constants of quantum electrodynamics (QED), \(\alpha \) [1,2,3], and quantum chromodynamics (QCD), \(\alpha _{\textrm{s}}\) [4, 5]; anomalous magnetic moment of electron, \(g_e\) [1, 6] and muon, \(g_\mu \) [7,8,9]; higher-order QED processes which have no classical analogy, such as recently measured light-by-light scattering [10, 11]. This paper explores a connection between QED and QCD virtual particles, which is predicted to give rise to a residual attractive force. This attractive force stems from gluon exchanges between strongly interacting degrees of freedom in off-shell photon.

The presence of strongly interacting degrees of freedom is well known to affect the value of \(\alpha \), \(g_e\) or \(g_\mu \) where it is quantified by the hadronic vacuum polarization (HVP) of photon propagator, \(\Pi ^{\textrm{had}}\). Baseline calculations of HVP were done by several authors [12,13,14]. Since the strong interaction turns non-perturbative due to the rise of \(\alpha _{\textrm{s}}\) at low energy, the low energy behavior of \(\Pi ^{\textrm{had}}\) is difficult to calculate using perturbative methods. Instead, it can be determined via the optical theorem, which connects \(\Pi ^{\textrm{had}}\) with the total cross section for the production of hadrons in \(e^+e^-\) collisions. Besides this approach, the \(\Pi ^{\textrm{had}}\) can be also determined from lattice QCD with increasing precision [9, 15].

The contributions to \(\Pi ^{\textrm{had}}\) at the momentum transferred (\(q^2\)) at the scale of hadron masses can be thought of as being due to exchanges of hadrons (predominantly mesons) in the loops [16]. For the \(q^2\) scale significantly below the hadron masses (\(q^2 \ll m_{\textrm{had}}\)), the relevant degrees of freedom in HVP are quarks and gluons. These are the same degrees of freedom as those responsible for the confinement of quarks in hadrons. One may, therefore, expect a negligible but nonzero confining force to act between the quarks in loops inside photons surrounding charged particles as well. Since photon has no direct coupling to gluon, it is the quark loops which are the basic strongly interacting degrees of freedom within photon. It is multi-gluon exchanges between quarks rather than exchanges within quark loops which seem to be relevant for the confinement [17]. Thus, multi-gluon exchanges between two quark loops should represent the basic contribution forming this force. An example of Feynman diagram contributing to this force is shown in Figure 1a. This ‘new’ force, which is, however, purely a Standard Model force, may be expected to have experimentally observable consequences and may obscure explanations of important observations at the scale of the gravity interaction, such as signatures of the dark matter.

In Sect. 2 of this paper, we calculate the asymptotic behavior of inter-quark-loop interactions. The contributions of these interactions to the scattering amplitude are then resummed in Sect. 3. The sign and strength of the interaction is discussed in Sect. 4. Observable consequences of this residual force are then discussed in Sect. 5.

Fig. 1
figure 1

Example of Feynman diagram for multiple gluon exchanges between two quark loops (a). Inter-quark-loop exchanges for 3-point quark loops without self-interactions (b, c), with self-interactions (d, e), and exchanges between 3- and 4-point quark loops (f). Only diagrams corresponding to non-abelian cases are shown (for the definition, see the text)

Full size image

2 Asymptotic behavior of inter-quark-loop interactions

The quark degrees of freedom inside photons surrounding two external charges may interact via the exchange of vector particles (photons, electroweak bosons, and gluons) and scalar particle (Higgs). Both electroweak boson and Higgs are too heavy to appear as effective degrees of freedom in the \(q^2 \ll m_{\textrm{had}}^2\) region. We will, therefore, discuss only the exchanges of photons between two quark loops (referred to as abelian case), and exchanges of gluons between two quark loops (referred to as non-abelian case). The latter is then the case where attractive confining forces may take place. The amplitude for these interactions can be expressed as

$$\begin{aligned} i \mathcal {M}^{\alpha \beta } &= (i\sqrt{a})^2 \frac{-i g^{\alpha \mu }}{q^2 +i\epsilon } Y_{\mu \nu }(q) \frac{-ig^{\nu \beta }}{q^2 + i\epsilon }, \end{aligned}$$
(1)
$$\begin{aligned} Y_{\mu \nu }(q)= & {} \int \prod _{i} \textrm{d}u_i A_{\mu ,\mu _1 \cdots \mu _{n-1}}(u_i,q) B^{\mu _1 \cdots \mu _{n-1},\nu _1 \cdots \nu _{n'-1}}(u_i, q) A'_{\nu ,\nu _1 \cdots \nu _{n'-1}}(u_i,q), \end{aligned}$$
(2)

where A and \(A'\) represent n- and \(n'\)-point quark loop amplitudes, respectively, B represents exchanges of boson fields, and \(u_i\) are unconstrained four-momenta (unconstrained four-momenta in quark loops are not explicitly written). We start the discussion with only the inter-loop exchanges. The intra-loop exchanges are mentioned later. Here, we also ignore the non-contracted wave function of incoming and outgoing charges, represented by bispinors and corresponding creation and annihilation operators. In the non-relativistic limit, these have an impact on the sign of the interaction, which is discussed in Sect. 4.

For the abelian case, \(n=n'\), and there are \(n-2\) momenta \(u_i\) connecting quark loops forming together multi-loop (n-loop) integral. For both abelian and non-abelian case the amplitude consists of at least three loops. The amplitudes A and \(A'\) also depend on quark masses, m, and B depends on intermediate boson mass, \(M\), which is also not explicitly written in (2). In the abelian case, \(M\) represents an infrared cutoff, which should be sent to zero at the end of the calculation. In the non-abelian case at low energies, \(M\) is the effective gluon mass stemming from non-perturbative aspects of infrared finite gluon propagator and should not be neglected (see, e.g., one of the latest calculations in Ref. [18]).

One loop amplitudes may be calculated using extensions of the original Passarino-Veltman (PV) reduction [19] or several new methods targeting mainly the LHC physics (see, e.g., review [20]). Calculations of multi-loop integrals remain an open problem though. Moreover, in the low-\(q^2\) regime, amplitudes such as the one in Fig. 1a are perturbative in \(\alpha \), but non-perturbative in the non-abelian part, meaning that adding more gluon exchanges may bring a larger contribution to the full amplitude. For those reasons, we will not attempt to provide a full solution of amplitudes at a given order, but we will discuss the general behavior of these amplitudes. In particular, we will explore the behavior of these amplitudes in the non-relativistic limit of \(q^2 \rightarrow 0\). For some of the calculations presented in the remainder of this section, we used software for automated calculation of loop integrals [21] interfaced with Python scripts allowing to sort or isolate various contributions in an efficient way.

The first possible contribution with \(n=n'=2\) and single gluon exchange between quark loops not captured by the formula (2) is forbidden due to color-charge conservation.

The first contribution we will analyze is therefore the contribution with \(n=n'=3\), that is the exchange between two 3-point quark loops, which appears both in the abelian and non-abelian case. There are two diagrams contributing, Fig. 1b and c, forming three-loop amplitudes \(Y_{\mu \nu }^b\) and \(Y_{\mu \nu }^c\), respectively. The space-time part of amplitude \(Y_{\mu \nu }^b\) is given by

$$\begin{aligned} Y_{\mu \nu }^b = \alpha _{(\textrm{s})}^2 \int \textrm{d}^4 u_1 & {} \int \textrm{d}^4 p \hspace{2.77771pt}\textrm{Tr} \bigg \{ \gamma _\mu \frac{i}{{p}\!\!\!/+{q}\!\!\!/-m} \gamma _{\mu _1} \frac{i}{{p}\!\!\!/+{u}\!\!\!/_1+{q}\!\!\!/-m} \gamma _{\mu _2} \frac{i}{{p}\!\!\!/ -m} \bigg \} \frac{-i g_{\mu _1 \nu _1}}{(u_1+q)^2-M^2} \frac{-i g_{\mu _2 \nu _2}}{u_1^2-M^2} \\{} & {} \int \textrm{d}^4 p' \hspace{2.77771pt}\textrm{Tr}\bigg \{\gamma _\nu \frac{i}{{p'}\!\!\!\!/ -m} \gamma _{\nu _2} \frac{i}{{p'}\!\!\!\!/+{u}\!\!\!/_1+{q}\!\!\!/-m} \gamma _{\nu _1} \frac{i}{{p'}\!\!\!\!/+{q}\!\!\!/-m} \bigg \}, \end{aligned}$$

where \(\alpha _{(\textrm{s})}\) is \(\alpha _{\textrm{s}}\) and \(\alpha \) for abelian and non-abelian case, respectively. The amplitudes \(Y_{\mu \nu }^b\) and \(Y_{\mu \nu }^c\) come out zero at various levels. First, they are zero due to Furry’s theorem (the odd number of bosons attached to the quark loop) for the abelian case as well as for the non-abelian case where the color structure does not protect the amplitude (trace over two Gell-Mann matrices is symmetrical under the exchange of color indices). Aside from that, the individual amplitudes may be calculated exactly in the limit \(q^2 \rightarrow 0\). When summed up, the two contributions come out with equal magnitudes but opposite signs, and they sum up to zero even at this level. Finally, each amplitude is equal to zero separately for \(M^2 \rightarrow 0\). For the non-abelian case, there are two more amplitudes with \(n=n'=3\), namely Fig. 1d, labeled \(Y_{\mu \nu }^d\), and Fig. 1e, labeled \(Y_{\mu \nu }^e\) (and amplitude with crossings). Again, these amplitudes come out zero due to Furry’s theorem. Furry’s theorem renders zero also amplitudes with \(n=3\) and \(n'>3\), such as \(Y_{\mu \nu }^f\) shown in Fig. 1f.

To further explore the basic behavior of non-Abelian amplitudes in \(q^2 \rightarrow 0\) limit, we step aside from Furry’s theorem and attempt to calculate individual amplitude \(Y_{\mu \nu }^f\). After evaluating two quark loops and performing contraction due to bosonic fields, one is left with approximately three thousand terms where a scalar product \(u_1\cdot u_2\) occurs inside PV coefficients. This makes it impossible to calculate the integral over \(u_i\) exactly using tools for automatic calculations. To obtain at least a partial result, we may isolate terms free of scalar products and perform the integration over \(u_i\) for these terms. This part of the result can be written in a form

$$\begin{aligned} Y_{\mu \nu }^f = g_{\mu \nu } \alpha _{\textrm{s}}^{\frac{5}{2}}\bigg [ \sum _{j=0}^{4} c_{0j} \log ^j\frac{\mu ^2}{M^2} + \frac{1}{\epsilon } \sum _{j=0}^{3} c_{1j} \log ^j\frac{\mu ^2}{M^2} + \frac{1}{\epsilon ^2}\sum _{j=0}^{2} c_{2j} \log ^j\frac{\mu ^2}{M^2} + \mathcal {O}(\epsilon ) \bigg ], \end{aligned}$$
(3)

where \(\epsilon \) and \(\mu ^2\) regularize UV divergences and \(c_{ij}=c_{ij}(m^2,M^2,\log (m^2/M^2),\textrm{PV}(m,M))\) are rational functions in \(m^2\) and \(M^2\), polynomial in \(\log (m^2/M^2)\) (up to the second order) and polynomial in PV functions (up to the second order). The maximal power of four of \(\log ({\mu ^2}/{M^2})\) is connected with the presence of four off-shell four-momenta that are integrated over. The maximal power of two of \(\log ({m^2}/{M^2})\) is connected with integrating over two four-momenta of intermediate bosons with effective mass \(M\). When performing the limit \(M\rightarrow 0\) on this rather complex result, one obtains zero. This infrared finite behavior is not directly visible from complex coefficients \(c_{ij}\) or from (3). In the following paragraphs, we will derive and generalize this infrared behavior for the analytical part of any two fermionic loop amplitude with \(n=n'\) and n even in the \(q^2 \rightarrow 0\) limit where \(n-1\) bosons are exchanged between the two loops. An example of such amplitude detailing the inter-quark-loop exchanges is shown in Fig. 2. For the abelian case, Furry’s theorem leaves nonzero only amplitudes with the odd number of boson exchanges between the two quark loops (n even, that is, even number of vertices with bosons when including also photon connecting the loop with external charge). For the non-abelian case, Furry’s theorem cannot be applied due to the color part of the amplitude. Consequently, we are left also with amplitudes having an even number of boson exchanges between the two quark loops, i.e., n odd. We perform the derivation only for n even.

Fig. 2
figure 2

Example diagram for n-1 boson exchanges between two n-point quark loops with boson four-momenta \(u_i\) and quark loop four-momenta p and \(p'\). Only a part of quark loops is shown here – the full quark loops, photons connecting loops with incoming and outgoing asymptotic states, as well as the asymptotic states as shown in Fig. 1 are suppressed here for legibility

Full size image

The Lorentz structure of one n-point fermionic loop calculation can be written as

$$\begin{aligned} A_{\mu _0,\mu _1 \cdots \mu _n}= & {} \int \textrm{Tr}\bigg \{ \prod _{i=0}^{n-1} \gamma _{\mu _i} \frac{i}{{p}\!\!\!/ + {u}\!\!\!/_i - m} \bigg \} \textrm{d}^4 p \nonumber \\{} & {} \sim \prod _{j=0}^{n-1} \sum _{i=0}^{n-1} u_{i\mu _j} + \sum _{\sigma (\mu )} \sum _{l=0}^{n/2-1} \bigg [ \prod _{l'=1}^{2l+\delta }\bigg ( \sum _{i=0}^{n-1} u_{i\mu _{n-l'}} \bigg ) \prod _{l'=1}^{(n-2l)/2} \bigg ( g_{\mu _{2l'-1}\mu _{2l'-2}} \bigg ) \bigg ], \end{aligned}$$
(4)

where \(u_n=0\), \(u_0=0\) for the case of \(q^2 \rightarrow 0\) limit, \(\delta =0\) and 1 for even and odd n, respectively, and sum with \(\sigma (\mu )\) encodes all possible permutations of Lorentz indices \(\mu _l\) (after expanding sums and products in (4), \(\mu _0\) should be relabeled to \(\mu \) to match the convention used in (2)). Equation (4) says that after performing the integration over p, the result is a sum of all possible combinations of vectors \(u_{i\mu _{l}}\) and metric tensors preserving the tensor structure of the original integral. Equation (4), therefore, captures only the Lorentz structure of the solution of one loop integral. The full solution is given by not written special functions multiplying each term in the sum (4) (e.g., PV functions or hypergeometric functions [22, 23]). Some terms in the full solution may contain scalar products \(u_i \cdot u_j\) as arguments of special functions. Some of these are non-analytical in \(u_i\cdot u_j\), other are analytical, can be expanded in a series, and contribute the terms of the same form as in r.h.s. of (4). We will further discuss only the analytical part of the one-loop solution.

After contracting the solution of quark loops with numerators of boson propagators, we are left with the numerator, \(\mathcal {N}\), of the multi-loop integral over \(u_i\), which has the following form

$$\begin{aligned} \sum _j \big [ g_{\mu \nu } c_1(u_i \cdot u_j, m) + u_{i\mu } u_{j\nu } c_2(u_i \cdot u_j, m) \big ], \end{aligned}$$
(5)

where j includes i and \(c_{1,2}\) are polynomials in scalar product \(u_i \cdot u_j\) and polynomials in PV functions with arguments \(u_i \cdot u_j, m\). When later integrating over a given four-momentum \(u_i\), this numerator can be written as a sum of terms proportional to

$$\begin{aligned} 1, ~ u_{i\alpha }, ~ u_{i\alpha }u_{i\beta }, ~ u_{i\alpha }u_{i\beta }u_{i\gamma }, ... , \end{aligned}$$
(6)

if excluding contributions non-analytical in \(u_i \cdot u_j\) as discussed earlier.

In the case of an odd number of boson exchanges between two quark loops, there is inevitably at least one factor with the following integral over \(u_i\),

$$\begin{aligned} \int \frac{\mathcal {N}}{ [(u_i+v)^2 - M^2][u_i^2 - M^2] } \textrm{d}^4 u_i, \end{aligned}$$
(7)

where \(v = \sum _{j,j \ne i} u_j\). This integration is then followed by at least one integral of the form

$$\begin{aligned} \int \frac{\mathcal {N}}{ u_j^2 - M^2 } \textrm{d}^4 u_j, ~ ~ \text {where } j \ne i. \end{aligned}$$
(8)

Performing the integral over \(u_i\), implies integrating (7) along with numerators (6). The result may be written for a given numerator \(\mathcal {N}\) as

$$\begin{aligned}{} & {} \mathcal {N} \sim 1 : \sim 1 + \frac{1}{\epsilon } + \log \frac{\mu ^2}{M^2}, \nonumber \\{} & {} \mathcal {N} \sim u_{i\alpha } : \sim v_\alpha \bigg ( 1 + \frac{1}{\epsilon } + \log \frac{\mu ^2}{M^2} \bigg ), \nonumber \\{} & {} \mathcal {N} \sim u_{i\alpha } u_{i\beta } : \sim \bigg ( v_\alpha v_\beta + M^2 g_{\alpha \beta } \bigg ) \bigg ( 1 + \frac{1}{\epsilon } + \log \frac{\mu ^2}{M^2} \bigg ), \nonumber \\{} & {} \mathcal {N} \sim u_{i\alpha } u_{i\beta } u_{i\gamma } : \sim \bigg ( v_\alpha v_\beta v_\gamma + \sum _{\sigma (\alpha \beta \gamma )} M^2 v_\alpha g_{\beta \gamma } \bigg ) \bigg ( 1 + \frac{1}{\epsilon } + \log \frac{\mu ^2}{M^2} \bigg ), \cdots . \end{aligned}$$
(9)

These are then integrated along with the denominator in (8) over each \(u_j\). Terms in the numerator with odd powers of \(u_j\) are equal to zero, and other terms lead to:

$$\begin{aligned}{} & {} \mathcal {N} \sim 1 : \sim M^2 \bigg ( 1 + \frac{1}{\epsilon } + \log \frac{\mu ^2}{M^2} \bigg ), \nonumber \\{} & {} \mathcal {N} \sim u_{j\alpha } u_{j\beta } : \sim M^4 g_{\alpha \beta } \bigg ( 1 + \frac{1}{\epsilon } + \log \frac{\mu ^2}{M^2} \bigg ), \nonumber \\{} & {} \mathcal {N} \sim u_{j\alpha } u_{j\beta } u_{j\gamma } u_{j\delta } : \sim M^6 \sum _{\sigma (\alpha \beta \gamma \delta )} g_{\alpha \beta } g_{\gamma \delta } \bigg ( 1 + \frac{1}{\epsilon } + \log \frac{\mu ^2}{M^2} \bigg ),\cdots , \end{aligned}$$
(10)

which may be generalized to

$$\begin{aligned} \mathcal {N} \sim \prod _{l=1}^{n} u_{j\mu _l} ~ ~ : ~ ~\sim M^2 \sum _{l=1}^{n/2} M^{n-2l} \sum _{\sigma (\mu )} \prod _{l'=1}^{(n-2l)/2}g_{\mu _{2l'-1}\mu _{2l'}} \bigg ( 1 + \frac{1}{\epsilon } + \log \frac{\mu ^2}{M^2} \bigg ). \end{aligned}$$
(11)

The full numerator on the l.h.s. of (10) and (11) also includes \(1/\epsilon \) and \(\log (\mu ^2/M^2)\) terms from (9). We will further concentrate only on the infrared part of the amplitude, which we label \(f_M\). It can then be written as

$$\begin{aligned} f_M(\alpha _s,M) \sim \alpha _s^{n-1} \big [ M^2( 1 + \log M^2 + \log ^2 M^2) + \mathcal {O}(M^4) \big ]. \end{aligned}$$
(12)

Since \(\lim _{M\rightarrow 0} (M^2 \log ^k(\mu ^2/M^2)) = 0\) for any finite k, the infrared part of the amplitude is finite and nonzero for the non-abelian case.

Apparently, the above shown solution describes only a certain sub-class of amplitudes. As said before, we left out non-abelian amplitudes with n odd. Further, intra-loop exchanges eliminate the applicability of Furry’s theorem even for the abelian amplitudes and may leave various other contributions. The additional contributions may also come from amplitudes with self-interactions. Nevertheless, this particular solution allows us to write down a general form of the solution for a sum of all configurations with bosonic interactions between two quark loops as

$$\begin{aligned} Y_{\mu \nu } = \bigg ( g_{\mu \nu } - \frac{q_\mu q_\nu }{q^2} \bigg ) \bigg (f_{\mathbb {A}}(\alpha ,m,q) + f_{\mathbb {N}}(\alpha ,m,q) \bigg )...\;\text {for abelian case,} \end{aligned}$$
(13)
$$\begin{aligned} Y_{\mu \nu } = \bigg ( g_{\mu \nu } - \frac{q_\mu q_\nu }{q^2} \bigg ) \bigg (f_M(\alpha _{\textrm{s}},M) + f_{\mathbb {A}}(\alpha _{\textrm{s}},m,M,q) + f_{\mathbb {N}}(\alpha _{\textrm{s}},m,M,q) \bigg )...\;\text {for non-abelian case,} \end{aligned}$$
(14)

where \(f_{\mathbb {A}}\) represents the function analytical in q, \(f_{\mathbb {N}}\) represents the non-analytical function (unsolved part of the amplitude), and \(f_M\) is finite and nonzero if \(M\) is nonzero. One may see this as a trivial solution especially since no further information about \(f_{\mathbb {A}}\) and \(f_{\mathbb {N}}\) is provided, but importantly, this solution encodes the finding that the nonzero M is a sufficient condition for a nonzero final amplitude in the non-relativistic limit (provided that \(f_{\mathbb {A}} \ne - f_M\), which is discussed in the next section). This condition is not met in the abelian case. As we only concentrate on the infrared part of the amplitude, the solution in (13) and (14) omits the short-distance ultraviolet part, which has a standard interpretation in terms of renormalizing bare quantities. While not providing much concrete information, the solution in (13) and (14) allows to resum the inter-quark-loop interactions, which is discussed in the next section.

It is also important to discuss the value of \(\alpha _{\textrm{s}}\) present in the amplitude. The \(\alpha _{\textrm{s}}\) has a Landau pole at \(q^2 = \Lambda _{\textrm{QCD}}^2\) at the leading order. However, the realistic estimate of \(\alpha _{\textrm{s}}\) for \(q^2 \rightarrow 0\) is \(\alpha _{\textrm{s}}(0) \approx 0.7-3\) [17]. Therefore, the value of \(\alpha _{\textrm{s}}(0)\) should not lead to a change in the presented conclusions.

3 Resumming inter-quark-loop interactions

The result presented in the previous section represents the calculation of the interaction between only two quark loops. To obtain a full result, we need to resum at the operator level all the contributions of this type, which may take place anywhere in the photon. As we could see, these contributions contain finite terms in the non-abelian case. These contributions describe the strong interaction between hadronic degrees of freedom, which may lead to a physical force. This physical force based on the strong interaction is then different than the force based on the electromagnetic interaction. For this reason, we should not resum these contributions together with the free propagator of electromagnetic interaction, but we should resum them separately.

If we label free photon propagator times \((-a)\) as 1/X and Y is defined in (2), we perform the following resummation,

$$\begin{aligned} \frac{1}{X} Y \frac{1}{X} +\frac{1}{X} Y \frac{1}{X} Y \frac{1}{X} + ... = \frac{1}{X} Y \bigg [ \frac{1}{X-Y}\bigg ]. \end{aligned}$$
(15)

Now, we can plug in the result at the eigenvalue level (14) and expand around \(q^2 = 0\) to perform the \(q^2 \rightarrow 0\) limit. When plugging (14) into (15) and expanding, one obtains

$$\begin{aligned} i\mathcal {M}^{\alpha \beta } = - g^{\alpha \beta }\bigg ( \frac{ia}{q^2} - \frac{1}{f_M(\alpha _{\textrm{s}},M) + f_{\mathbb {A},0}(\alpha _{\textrm{s}},m,M) + f_{\mathbb {N}}(\alpha _{\textrm{s}},m,M,q)} + \mathcal {O}(q) \bigg ), \end{aligned}$$
(16)

where \(f_{\mathbb {A},0}(\alpha _{\textrm{s}},m,M)\) is the zeroth term in the expansion of \(f_{\mathbb {A}}(q)\). The \(1/q^2\) behavior of the amplitude implies 1/r behavior of the resulting non-relativistic potential. It is noticeable that the non-analytical part, which is undetermined and which may be given, e.g., by a sum of polylogarithms, stands in the denominator along with \(f_M\). If \(M\) is nonzero, any behavior of \(f_{\mathbb {N}}\) with \(q^2 \rightarrow 0\) will keep unchanged the \(1/q^2\) behavior of the amplitude. Namely, for \(f_{\mathbb {N}} \rightarrow \infty \) in the limit of \(q^2 \rightarrow 0\), all the terms but the first one vanish. For \(f_{\mathbb {N}}\) finite or zero in the limit of \(q^2 \rightarrow 0\), only the first term in (16) is divergent, and the second finite term represents a finite contribution leading to a local correction of the 1/r behavior of the potential. Terms with \(1/\sqrt{q^2}\) or \(\ln q^2\), which are divergent for \(q^2 \rightarrow 0\) and which would give rise to corrections of potential behaving as \(1/r^2\) or \(1/r^3\), respectively [24], are absent.

We should note that \(f_{\mathbb {N}}\) is present in the numerator not only with \(f_M\) but also with \(f_{\mathbb {A},0}\). We did not prove in general that \(f_{\mathbb {A},0} \ne -f_M\), which would lead to a cancellation and possible deviation from 1/r behavior due to the presence of \(1/f_{\mathbb {N}}(q)\) in the amplitude, but we verified that \(f_{\mathbb {A},0} \ne -f_M\) for amplitudes discussed in Sect. 2.

There are two other observations that can be made. First, one can compare the amplitude (16) to the leading order scattering amplitude for the photon exchange between two fermions of the same electric charge, which reads \(i\mathcal {M}^{\alpha \beta } = g^{\alpha \beta }\frac{ia}{q^2}\) (still ignoring the description of the observable initial and final state which will be discussed in the next section). From this comparison, one can see that there is a change in the sign between the two amplitudes which originates in the resummation (15) and subsequent expansion in q. The details about the sign and strength of the interaction are then discussed in the next section. Second, one can easily verify that the sign of Y influences only the sign of the second and other terms in (16). The sign of contributions to Y is indeed alternating due to the presence of different number of interacting fields, due to the color algebra, or due to the presence of amplitudes with crossings.

4 Sign and strength of the interaction

We will now turn attention to the problem of the sign of the interaction. First, we will briefly recap the origin of the sign in two basic theories – in Yukawa and electromagnetic theory. As this should be well known from classical textbooks (see, e.g., [25]), we will recap that in the form of a table where the multiplicative sign factors and their source are shown for interaction between two fermions. The origin of sign of interaction for fermion-fermion (\(ff\)), fermion-anti-fermion (\(f\bar{f}\)), and two different fermions (\(f_{1}^{-}f_{2}^{+}\) and \(f_{1}^{-}\bar{f}_{2}^{-}\)) is summarized in Table 1.

Table 1 The origin of sign in Yukawa and electromagnetic interaction for fermion-fermion (\(ff\)), fermion-anti-fermion (\(f\bar{f}\)), and two different fermions (\(f_{1}^{-}f_{2}^{+}\) and \(f_{1}^{-}\bar{f}_{2}^{-}\)). Details are provided in the text
Full size table

The fermion current sign comes from the normalization of bispinors in \(\bar{\Psi }\Psi \) and \(\bar{\Psi }\gamma ^\mu \Psi \) currents in Yukawa and electromagnetic interaction, respectively. In the non-relativistic limit, only the term with \(\gamma ^0\) remains due to the definition of bispinors, which compensates the sign from the normalization of bispinor in the case of electromagnetic interaction between fermion-antifermion. Since only \(\gamma ^0\) term is active in fermion currents, only \(g_{00}\) term remains in the photon propagator, leading to the repulsive instead of attractive interaction in the case of electromagnetic interaction between \(ff\) pair (intermediate particle sign). While the minus sign from the normal product of creation and annihilation operators compensates the sign from bispinor normalization in the Yukawa case, it leads to the total attractive force in the case of \(f\bar{f}\) electromagnetic interaction. Since f and \(\bar{f}\) represent two states of one wave function, the sign of the charge in the coupling remains the same and it is the order of creation and annihilation operators in the normal product which dictates the alternating sign of the corresponding Coulomb potential. In the case of \(f_{1}^{-}f_{2}^{+}\), the order of creation and annihilation operators is the same as in the \(ff\) case. The order is the same even in the case of exchanging \(f_2\) by its antiparticle since \(\bar{f_2}\) is not an antiparticle to \(f_1\). In the case of \(f_{1}^{-}f_{2}^{+}\) and \(f_{1}^{-}\bar{f}_{2}^{-}\), the total sign of electromagnetic interaction is dictated by the sign of charge in the coupling. The factor \((-1)\) in front of the total sign comes from \(q^2\) being space-like, but in general, it is only the relative sign that can be determined while the total sign is a matter of convention.

In the case of residual strong interaction, the signs due to fermion current and normal product do not change with respect to the electromagnetic interaction. Neither changes the sign due to the presence of \(g_{00}\) in the non-relativistic limit. What changes is the overall sign, which is due to resummation as discussed in Section 3. Since the physical interaction is due to the strong interaction which does not distinguish charge and which probes the strongly interacting degrees of freedom in the wave function of incoming fermions, the sign of the electromagnetic charge is not relevant. The sign in coupling is thus the same as in the Yukawa case leading to the attractive interaction for all the configurations except for the case of \(f\bar{f}\) pair where it is repulsive. This repulsive behavior in \(f\bar{f}\) case has the same origin as the attractive behavior in the electromagnetic interaction of \(f\bar{f}\), which is due to the description of physical, asymptotic states as discussed in previous paragraphs. This prediction represents a distinguishable behavior of the residual force.

Since our interaction probes only hadronic degrees of freedom in the fermion’s wave function, the magnitude of coupling is given by the hadronic part of electromagnetic coupling. The interaction strength of strongly interacting part of electromagnetic fields can be determined from the leading order expansion of the electromagnetic coupling constant \(\alpha (q^2)\) which reads

$$\begin{aligned} \alpha (q^2) & = {} \frac{\alpha _0}{1 - \varDelta \alpha (q^2)} = \alpha _0 [ 1 + \varDelta \alpha (q^2)_{\textrm{had}} \nonumber \\{} & \quad {} + \varDelta \alpha (q^2)_{\textrm{lep}} + \varDelta \alpha (q^2)_{\textrm{top}} + \mathcal {O}(\varDelta \alpha ^2)]. \end{aligned}$$
(17)

This gives an estimate of the strength of interaction mediated by a photon containing quark loops, which is

$$\begin{aligned} \alpha _0 \varDelta \alpha (q^2)_{\textrm{had}} \equiv a(q^2). \end{aligned}$$
(18)

The \(a(q^2)\) is then the value of coupling to be used in resummed expressions (16).

5 Estimates of the magnitude of residual force at macroscopic scales

As discussed in previous sections, the residual force may be expected to act as an imperceptible, attractive component of the Coulomb potential, which does not distinguish the sign of the charged object. It may therefore be detectable at macroscopic scales as a difference in the magnitude of the repulsive Coulomb force between two charged objects with equal charges \(e_1=e_2=e\) and the magnitude of the attractive Coulomb force between two opposite charges with \(e_1=-e_2=e\).

We estimate the ratio of magnitudes of the residual force and electromagnetic force, \(a/\alpha _0\), using the values of \(\varDelta \alpha _{\textrm{had}}\) based on numerical calculations from [26], which use extrapolations from hadron measurements as briefly discussed in Section 4 and detailed in Ref. [27]. The result is \(10^{-20}\) and \(10^{-26}\) for the eV scale and meV scale, respectively. That is, for these scales, the ratio is not extremely far from the ratio of the magnitude of gravity and electromagnetic force for two protons (which is \(10^{-36}\)). We should stress that the presented numbers provide only a basic estimate based on extrapolations of trends in the evolution of \(\varDelta \alpha _{\textrm{had}}\) done over several orders of magnitude, and as such, they should be taken with a reserve. Precise estimates of \(\varDelta \alpha _{\textrm{had}}\) based on lattice QCD would be needed to build a firm quantitative statement. Nevertheless, we can compare these estimates with the precision of current terrestrial experiments.

Assuming that terrestrial experiments with macroscopic electromagnetic forces can be done around the energy level of eV, detecting the residual force at terrestrial conditions would require precision of the measurement of electromagnetic processes at the level of \(10^{-20}\). The relative uncertainty of the knowledge of elementary charge and permittivity of vacuum is at the level of \(10^{-10}\) [28]. The precision of the measurement of Coulomb force for elementary charges is at a similar level [29]. Thus, it seems impossible at the present time to verify or falsify the presence of residual force by terrestrial experiments.

While it is not possible to measure the residual force at terrestrial scales, it may play a role in large-scale astrophysical objects due to the presence of QED plasma in thermal equilibrium. Two examples of those objects are the warm ionized medium (WIM) in interstellar space [30,31,32] and intracluster medium (ICM) in clusters of galaxies [33, 34]. Those media are composed predominantly of ionized hydrogen. They are macroscopically electrically neutral due to Coulombic interactions between ions and electrons acting within a Debye volume defined by the Debye length, \(\lambda _{\textrm{D}}\). The screened electrostatic potential is given by [35]

$$\begin{aligned} \Phi _{\textrm{deb}}(r) = \frac{1}{4\pi \epsilon _0} \frac{e}{r} \exp \bigg ( -\frac{\sqrt{2}r}{\lambda _{\textrm{D}}} \bigg ), \end{aligned}$$
(19)

where e is the unit charge. While the QED interactions among positive and negative charges lead to net-zero Coulomb force among plasma constituents, the residual force is attractive for both the positive and negative charges. The net-zero Coulomb force may, therefore, allow the residual force to be macroscopically distinguishable.

The estimate of the magnitude of the residual force between a given H\(^+\) ion and other ions in the Debye volume can be obtained by replacing \(\alpha _0\) with \(a(q^2)\) in a classical calculation of force due to potential (19) which is

$$\begin{aligned} F_{\textrm{res}}(\lambda _{\textrm{D}})& = {} 2 n_i e \int ^{\lambda _{\textrm{D}}} \frac{\partial \Phi _{\textrm{deb}}(r)}{\partial r} \textrm{d}^3 x \bigg |_{\alpha _0 \rightarrow a(q^2)} \nonumber \\ &= {} 4 \pi [\sqrt{2} - (1 \! + \! \sqrt{2})\exp (-\sqrt{2})] n_i a(q^2) \lambda _{\textrm{D}} \hbar c, ~ ~ \end{aligned}$$
(20)

where \(n_i\) is the number density of the ionized medium. The magnitude of residual force may be compared with the magnitude of gravitational force within the same volume

$$\begin{aligned} F_{\textrm{G}}(\lambda _{\textrm{D}}) = 4 \pi n_i \frac{m^2_{\textrm{H}^+}}{m_{\textrm{Pl}}^2} \lambda _{\textrm{D}} \hbar c , \end{aligned}$$
(21)

where \(m_{\textrm{Pl}}\) is the Planck mass and \(m_{\textrm{H}^+}\) is the proton mass.

For WIM, where \(n_i \sim 10^{-1}\) cm\(^{-3}\), \(\lambda _{\textrm{D}} \sim 10\) m, and temperature \(T \sim 10^{4}\) K (\(q^2\) at the scale of eV\(^{2}\)), the magnitude of \(F_{\textrm{res}}\) is \(10^{-46}\) GeV/fm. For ICM, where \(n_i \sim 10^{-3}\) cm\(^{-3}\), \(\lambda _{\textrm{D}} \sim 10^{3}\) m, and \(T \sim 10^7 - 10^8\) K (\(q^2\) at the scale of keV\(^{2}\)), the order of magnitude of \(F_{\textrm{res}}\) is \(10^{-40}\) GeV/fm. The order of magnitude of \(F_{\textrm{G}}\) is \(10^{-62}\) and \(10^{-61}\) GeV/fm for WIM and ICM, respectively. While the values of \(F_{\textrm{res}}\) should be taken with reserve as they are based on extrapolations of \(\varDelta \alpha _{\textrm{had}}\), as previously discussed, they suggest that the magnitude of \(F_{\textrm{G}}\) might be significantly smaller, implying that the presence of attractive force may be stronger than the self-gravity within the Debye volume. This might then have consequences on the understanding of the dynamics of these media.Footnote 1

6 Discussion and conclusion

Residual forces of electromagnetism and strong interaction give rise to many structures in nature, but their description by a theory from first principles is not simple. Here, we proposed a new residual force, which was not discussed previously in the literature, the residual QCD force in electromagnetically interacting fields. We calculated the asymptotic behavior of relevant scattering amplitudes, performed their resummation, and analyzed the sign of the resulting interaction. The non-relativistic potential of the residual force exhibits 1/r behavior and it is attractive for all particles except for the interaction between fermion and anti-fermion. Then, we calculated the primary experimentally observable consequences of this Standard Model force. While imperceptible at the scale of terrestrial experiments, the presence of this force may have an impact at astrophysical scale, in particular, it may influence the dynamics of the warm ionized medium present in galaxies and the intergalactic space.

Better quantitative understanding of the behavior of hadronic vacuum polarization, e.g., based on lattice QCD calculations, may bring further insight. Including an attractive force of this kind into the simulation of large-scale structures in the universe that involve QED plasma may also be fruitful. The presented paper opens the door for these investigations.