Renormalisable Non-Local Quark–Gluon Interaction: Mass Gap, Chiral Symmetry Breaking and Scale Invariance

Abstract

We derive a Nambu–Jona-Lasinio (NJL) model from a non-local gauge theory and show that it has confining properties at low energies. In particular, we present an extended approach to non-local QCD and a complete revision of the technique of Bender, Milton and Savage applied to non-local theories, providing a set of Dyson–Schwinger equations in differential form. In the local case, we obtain closed-form solutions in the simplest case of the scalar field and extend it to the Yang–Mills field. In general, for non-local theories, we use a perturbative technique and a Fourier series and show how higher-order harmonics are heavily damped due to the presence of the non-local factor. The spectrum of the theory is analysed for the non-local Yang–Mills sector and found to be in agreement with the local results on the lattice in the limit of the non-locality mass parameter running to infinity. In the non-local case, we confine ourselves to a non-locality mass that is sufficiently large compared to the mass scale arising from the integration of the Dyson–Schwinger equations. Such a choice results in good agreement, in the proper limit, with the spectrum of the local theory. We derive a gap equation for the fermions in the theory that gives some indication of quark confinement in the non-local NJL case as well. Confinement seems to be a rather ubiquitous effect that removes some degrees of freedom in the original action, favouring the appearance of new observable states, as seen, e.g., for quantum chromodynamics at lower energies.

1. Introduction

A lesson one can take from string theory is that strings are non-local objects [1,2,3,4,5,6,7,8,9,10,11,12]. Indeed, this lesson is the original motivation for the weakly non-local field theory to devise a novel pathway for possible UV-regularised theories inspired by string field theory [1,13], as investigated in several Refs. since the 1990s [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39]. Later on, this pathway became an alternative to address the divergence and the hierarchy problems in the Standard Model (SM) via generalising the kinetic energy operators that are of second order in the derivative to an infinite series of higher-order derivatives, suppressed by the scale $M^2$ of non-locality [14,40,41]. Such modifications of the kinetic energy sector in introducing higher-order derivatives are free of ghosts [33] (the unitarity issues are well addressed due to a certain prescription [42,43,44,45,46]) and also cure the Higgs vacuum instability problem of the SM Higgs [47], as analysed by one of the authors. It was shown that the $\beta$ function reaches a conformal limit, resolving the Landau pole issue in Quantum Field Theory [48]. Therefore, by capturing infinite derivatives by an exponential of an entire function, we obtain a softened ultraviolet (UV) behaviour in the most suitable manner, without the cost of introducing any new degrees of freedom that contribute to the particle mass spectrum, since there are no new poles in the propagators of infinite-derivative extensions. (For the astrophysical implications, dimensional transmutation and dark matter and dark energy phenomenology in these theories, see Refs. [49,50,51]). A bound on the scale of non-locality from observations from LHC and dark matter physics is $M\ge O\left(10\right)\mathrm{TeV}$ [41,50]. Moreover, the non-local theory leads to interesting implications for proton decay and Grand Unified Theories [52], as well as for braneworld models [53]. In addition, strongly coupled non-perturbative regimes, exact $\beta$ functions and the conditions of confinement in higher-derivative non-local theories are being actively investigated [54,55,56,57,58,59]. The results obtained so far show that the effect of the non-locality in the strong coupling limits is the dilution of any mass gaps (that may be present in the theory) in the UV regime that the system generates. (In the context of gravity theories, one can get rid of classical singularities, such as black hole singularities [19,24,60,61,62,63,64,65,66,67,68,69] and cosmological singularities [70,71,72,73,74,75,76]. Recently, false vacuum tunnelling studies were carried out in Ref. [77]. Relations to string theory are found in Refs. [78,79,80] (for an overview, cf. Ref. [81])).

In order to render charge renormalisation finite, infinite-derivative terms should be introduced for the fermion fields as well [48]. An analogous mechanism is the well-known Pauli–Villars regularisation scheme, where a mass dependence $\Lambda$ is introduced as the cut-off. The infinite-derivative approach is to promote the Pauli–Villars cut-off $\Lambda$ to the non-local energy scale M of the theory. In the infinite-derivative case, the non-local energy scale M might play the role of the ultraviolet cut-off. Recently, the infinite-derivative model for QED and Yang–Mills has been reconsidered in view of its generalisation to the SM [47,48]. (See also Refs. [41,50] for a few phenomenological applications to LHC physics and dark matter physics). Indeed, this model leads to a theory that is naturally free of quadratic divergencies, thus providing an alternative way to the solution of the hierarchy problem [41]. Higher-dimensional operators containing new interactions naturally appear in higher-derivative theories with a non-Abelian gauge structure. These operators soften quantum corrections in the UV regime and extinguish divergencies in radiative corrections. Nevertheless, as can be easily understood by power-counting arguments, the new higher-dimensional operators do not break renormalisability [48]. This is due to the improved ultraviolet behaviour of the bosonic propagator $P\left(k\right)$ in the deep Euclidean region, which scales as $P\left(k\right)\sim exp(-k^2/M^2)/k^4$, instead of the usual propagator scaling in the local case as $P\left(k\right)\sim 1/k^2$ for $k^2\to \infty$ [48].

Note that the presence of the non-local energy scale M associated with the infinite-derivative term manifestly breaks (at the classical level) the conformal symmetry of the unbroken gauge sector. Therefore, one may wonder whether this term can also trigger (dynamically) chiral symmetry breaking at low energy, or in other words, whether the fermion field could dynamically obtain a mass m satisfying the condition $m<M$. The aim of the present paper is to investigate this issue by analysing a general class of renormalisable models containing infinite-derivative terms. In this paper, we will show that a non-vanishing mass term for the fermion field can indeed be generated, depending on the kind of interaction at hand, as a solution of the mass gap equation. The fermion mass can be predicted, and it turns out to be a function of the energy scale M. The effect of non-locality is to move a possible violation of micro-causality to the region beyond the non-locality mass scale M, making them possibly unobservable. On the other hand, diagrammatic techniques a la Feynman cannot properly work in this context, which makes our approach through solutions to the Dyson–Schwinger set of equations more appealing.

This paper is organised as follows. In Section 2, we introduce the infinite-derivative $SU\left(N\right)$ gauge theory that we aim to study. In Section 3, we derive a set of Dyson–Schwinger equations for non-local QCD. We solve these equations in a perturbative manner by noting that higher harmonics are heavily damped by the non-local factor. In Section 4, the spectrum of the theory is analysed for the non-local Yang–Mills sector and found to be in agreement with the local results on the lattice in the limit of the non-local scale running to infinity. In the non-local case, we confine ourselves to a non-local scale that is sufficiently large with respect to the mass scale arising from the integration of the Dyson–Schwinger equations. In Section 5, we derive the gap equation for the fermion in the theory and show that an identical argument to that given in Refs. [58,59] can be applied here, giving some indication of quark confinement in the non-local case as well. Section 6 contains our conclusions. In Appendix A, we derive the Dyson–Schwinger equations for a ${\varphi }^4$ theory in differential form. In Appendix B, we calculate the scalar two-point function.

2. Infinite-Derivative SU(N)

The introduction of infinite-derivative terms in the Lagrangian is based on an approach by Lee and Wick [82,83] (cf. also Ref. [84]). The Lee–Wick approach, understood as the first terms in a series expansion, leads to a non-local infinite-derivative approach to $SU\left(N\right)$ gauge theories. In infinite-derivative $SU\left(N\right)$ with massless fermions, we start with the Lagrangian [48]:

$${\mathcal{L}}^{\prime }=\overline{\psi }{\gamma }^{\mu }(i{\partial }_{\mu }-g_s{A}_{\mu }^a{T}_a)\psi -\frac{1}{4}{F}_{\mu \nu }^a{e}^{-f\left(D^2\right)}{F}_a^{\mu \nu }-\frac{1}{2\xi }\left({\partial }_{\mu }{A}_a^{\mu }\right)e^{-f\left(D^2\right)}\left({\partial }_{\nu }{A}_a^{\nu }\right)+j_a^{\mu }{A}_{\mu }^a$$

(1)

where we assume the non-locality to be only in the gauge sector and $F_{\mu \nu }^a={\partial }_{\mu }{A}_{\nu }^a-{\partial }_{\nu }{A}_{\mu }^a-g_s{f}_{abc}{A}_{\mu }^b{A}_{\nu }^c$. Both the Yang–Mills term and the gauge-fixing term are delocalised by the infinite-derivative exponential. A source term $j_a^{\mu }{A}_{\mu }^a$ is added in order to use the Lagrangian for the generating functional. As an example, we can use $f\left(D^2\right)=D^2/M^2$, where $D_{\mu }^{ab}={\delta }^{ab}{\partial }_{\mu }-i{g}_s{A}_{\mu }^c{\left(T_c\right)}^{ab}$ is the covariant derivative in the adjoint representation. For a large non-local scale $M^2$, it was shown in Refs. [54,55] that $e^{f\left(D^2\right)}$ can be approximated by $e^{f(\square )}$, where $\square ={\partial }^2$. Therefore, one can start with

$$\mathcal{L}=\overline{\psi }{\gamma }^{\mu }(i{\partial }_{\mu }-g_s{A}_{\mu }^a{T}_a)\psi -\frac{1}{4}{F}_{\mu \nu }^a{e}^{-f(\square )}{F}_a^{\mu \nu }-\frac{1}{2\xi }\left({\partial }_{\mu }{A}_a^{\mu }\right)e^{-f(\square )}\left({\partial }_{\nu }{A}_a^{\nu }\right)+j_a^{\mu }{A}_{\mu }^a.$$

(2)

Dealing with the non-locality, we use a redefinition of fields by employing

$$A_{\mu }^a=e^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^a$$

(3)

which leads to

$$\mathcal{L}=\overline{\psi }{\gamma }^{\mu }(i{\partial }_{\mu }-g_s\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^a\right]T_a)\psi +{\mathcal{L}}_{\mathrm{YM}+},$$

(4)

where

$${\mathcal{L}}_{\mathrm{YM}+}=-\frac{1}{4}{\underline{F}}_{\mu \nu }^a{\underline{F}}_a^{\mu \nu }-\frac{1}{2\xi }\left[{\partial }_{\mu }{\underline{A}}_a^{\mu }\right]\left[{\partial }_{\nu }{\underline{A}}_a^{\nu }\right]+j_a^{\mu }\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^a\right].$$

(5)

The underline stands for the field redefinition. Here and in the following, the square brackets restrict the operational range of the differential operators included. Using the fact that the Lagrangian is determined up to a total divergence, one can use integration by parts to obtain

$$F_{\mu \nu }^a{e}^{-f(\square )}{F}_a^{\mu \nu }=F_{\mu \nu }^a{e}^{-\frac{1}{2}f(\square )}{e}^{-\frac{1}{2}f(\square )}{F}_a^{\mu \nu }=\left[e^{-\frac{1}{2}f(\square )}{F}_{\mu \nu }^a\right]\left[e^{-\frac{1}{2}f(\square )}{F}_a^{\mu \nu }\right]=:{\underline{F}}_{\mu \nu }^a{\underline{F}}_a^{\mu \nu }$$

(6)

up to a total divergence (named $div$ in the following and vanishing in the action integral), where $f[\square g]=[\square f]g+div$, and (by induction) $f\left[{\square }^{n}g\right]=\left[{\square }^{n}f\right]g+div$ is used to share the non-locality equally between the two field strength tensors. Explicitly, one has

$${\underline{F}}_{\mu \nu }^a={\partial }_{\mu }{\underline{A}}_{\nu }^a-{\partial }_{\nu }{\underline{A}}_{\mu }^a+g_s{f}_{abc}{e}^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^b\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\nu }^c\right].$$

(7)

For the modified Yang–Mills and gauge-fixing Lagrangian, one obtains

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\mathrm{YM}+}& =& -\frac{1}{4}[{\partial }_{\mu }{\underline{A}}_{\nu }^a-{\partial }_{\nu }{\underline{A}}_{\mu }^a][{\partial }^{\mu }{\underline{A}}_a^{\nu }-{\partial }^{\nu }{\underline{A}}_a^{\mu }]-\frac{1}{2\xi }\left[{\partial }_{\mu }{\underline{A}}_a^{\mu }\right]\left[{\partial }_{\nu }{\underline{A}}_a^{\nu }\right]+j_a^{\mu }\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^a\right]\phantom{(}\hfill \\ & & \phantom{(}-\frac{1}{2}{g}_s{f}_{abc}[{\partial }_{\mu }{\underline{A}}_{\nu }^a-{\partial }_{\nu }{\underline{A}}_{\mu }^a]e^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_b^{\mu }\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_c^{\nu }\right]\phantom{(}\hfill \\ & & \phantom{(}-\frac{1}{4}{g}_s^2{f}_{abc}{f}_{ade}\left[e^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_b^{\mu }\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_c^{\nu }\right]\right]\left[e^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^d\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\nu }^e\right]\right]\hfill \\ & =& -\frac{1}{2}\left[{\partial }_{\mu }{\underline{A}}_{\nu }^a\right][{\partial }^{\mu }{\underline{A}}_a^{\nu }-{\partial }^{\nu }{\underline{A}}_a^{\mu }]-\frac{1}{2\xi }\left[{\partial }_{\mu }{\underline{A}}_a^{\mu }\right]\left[{\partial }_{\nu }{\underline{A}}_a^{\nu }\right]-i{j}_a^{\mu }\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^a\right]\phantom{(}\hfill \\ & & \phantom{(}-\frac{1}{2}{g}_s{f}_{abc}[{\partial }_{\mu }{\underline{A}}_{\nu }^a-{\partial }_{\nu }{\underline{A}}_{\mu }^a]e^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_b^{\mu }\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_c^{\nu }\right]\phantom{(}\hfill \\ & & \phantom{(}-\frac{1}{4}{g}_s^2{f}_{abc}{f}_{ade}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_b^{\mu }\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_c^{\nu }\right]e^{-f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^d\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\nu }^e\right]+div,\hfill \end{array}$$

(8)

where the first term has been simplified, and integrations by parts have been used to move the exponential derivatives to the central position. This can also be reverted if necessary, i.e., if the derivatives are inappropriate for the variation in the fields. An ordinary integration by parts in the first line for the second part of $\mathcal{S}=\int \mathcal{L}{d}^{4}x={\mathcal{S}}_f+{\mathcal{S}}_{\mathrm{YM}+}$ gives

$$\begin{array}{ccc}\hfill {\mathcal{S}}_{\mathrm{YM}+}& =& \int d^{4}x[\frac{1}{2}{\underline{A}}_{\mu }^a\left(x\right)(\square {\eta }^{\mu \nu }-{\partial }^{\mu }{\partial }^{\nu }){\underline{A}}_{\nu }^a\left(x\right)+\frac{1}{2\xi }{\underline{A}}_{\mu }^a\left(x\right){\partial }^{\mu }{\partial }^{\nu }{\underline{A}}_{\nu }^a\left(x\right)+j_a^{\mu }\left(x\right)\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^a\left(x\right)\right]\phantom{(}\hfill \\ & & \phantom{(}-\frac{1}{2}{g}_s{f}_{abc}[{\partial }_{\mu }{\underline{A}}_{\nu }^a\left(x\right)-{\partial }_{\nu }{\underline{A}}_{\mu }^a\left(x\right)]e^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_b^{\mu }\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_c^{\nu }\left(x\right)\right]\phantom{(}\hfill \\ & & \phantom{(}-\frac{1}{4}{g}_s^2{f}_{abc}{f}_{cde}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_a^{\mu }\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_b^{\nu }\left(x\right)\right]e^{-f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^d\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\nu }^e\left(x\right)\right]].\hfill \end{array}$$

(9)

The Euler–Lagrange equation $\delta \mathcal{S}/\delta {\underline{A}}_{\mu }^a=0$ for the Yang–Mills field reads

$$\begin{array}{c}{\displaystyle g_s{e}^{\frac{1}{2}f(\square )}\overline{\psi }{\gamma }^{\mu }{T}_a\psi }\hfill \\ =(\square {\eta }^{\mu \nu }-{\partial }^{\mu }{\partial }^{\nu }){\underline{A}}_{\nu }^a+\frac{1}{\xi }{\partial }^{\mu }{\partial }^{\nu }{\underline{A}}_{\nu }^a+e^{\frac{1}{2}f(\square )}{j}_a^{\mu }-g_s{f}_{abc}{e}^{-\frac{1}{2}f(\square )}{\partial }_{\nu }\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_b^{\mu }\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_c^{\nu }\right]\phantom{(}\hfill \\ \phantom{(}-g_s{f}_{abc}{e}^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\nu }^b\right]e^{-\frac{1}{2}f(\square )}({\partial }^{\mu }{\underline{A}}_c^{\nu }-{\partial }^{\nu }{\underline{A}}_c^{\mu })+\phantom{(}\hfill \\ \phantom{(}-g_s^2{f}_{abc}{f}_{cde}{e}^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_{\nu }^b\right]e^{-f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_d^{\mu }\right]\left[e^{\frac{1}{2}f(\square )}{\underline{A}}_e^{\nu }\right].\hfill \end{array}$$

(10)

In the following, we use the Feynman gauge $\xi =1$ to further simplify the first line. Applying the formalism of Bender, Milton and Savage [85] to the Euler–Lagrange Equation (10) leads to the tower of Dyson–Schwinger equations, which is formulated in terms of multiple-point Green functions. This approach is displayed in Appendix A and Appendix B for a local ${\varphi }^4$ theory.

3. Solution of Dyson–Schwinger Equations

The creation of the tower of Dyson–Schwinger equations is explained in detail in Appendix B of Ref. [54]. We do not repeat it here to avoid reprinting previously published material. Instead, we give some guidelines. To start with, one takes the expectation value of the Euler–Lagrange Equation (10), weighted by the generating functional

$$Z\left[j\right]=\int \left[d\underline{A}\right]exp\left(i\int {\mathcal{L}}_{\mathrm{YM}+}{d}^{4}x\right),$$

(11)

where the current j is already contained in ${\mathcal{L}}_{\mathrm{YM}+}$, indicated by the plus sign in the index. One has $\langle e^{\frac{1}{2}f\left({\square }_x\right)}{\underline{A}}_{\mu }^a\left(x\right)\rangle =:Z\left[j\right]e^{\frac{1}{2}f\left({\square }_x\right)}{\underline{G}}_{1\mu }^{\left(j\right)a}\left(x\right)$. By either applying an additional partial derivative like $\langle e^{\frac{1}{2}f\left({\square }_x\right)}{\partial }_{\nu }^x{\underline{A}}_{\mu }^a\left(x\right)\rangle =Z\left[j\right]e^{\frac{1}{2}f\left({\square }_x\right)}{\partial }_{\nu }^x{\underline{G}}_{1\mu }^{\left(j\right)a}\left(x\right)$ or calculating a variation with respect to, e.g., $j_b^{\nu }\left(x^{\prime }\right)$,

$$\begin{array}{c}{\displaystyle \langle e^{\frac{1}{2}f\left({\square }_x\right)}{\underline{A}}_{\mu }^a\left(x\right)e^{\frac{1}{2}f\left({\square }_{x^{\prime }}\right)}{\underline{A}}_{\nu }^b\left(x^{\prime }\right)\rangle }\hfill \\ \hfill =Z\left[j\right]e^{\frac{1}{2}f\left({\square }_x\right)}{\underline{G}}_{2\mu \nu }^{\left(j\right)ab}(x,x^{\prime })+Z\left[j\right]e^{\frac{1}{2}f\left({\square }_x\right)}{\underline{G}}_{1\mu }^{\left(j\right)a}\left(x\right)e^{\frac{1}{2}f\left({\square }_{x^{\prime }}\right)}{\underline{G}}_{1\nu }^{\left(j\right)b}\left(x^{\prime }\right),\end{array}$$

(12)

one finally obtains

$$\begin{array}{c}{\displaystyle \square {\underline{G}}_{1\mu }^{\left(j\right)a}\left(x\right)+e^{\frac{1}{2}f(\square )}{j}_{\mu }^a\left(x\right)=}\hfill \\ =g_s{f}_{abc}\{e^{-\frac{1}{2}f(\square )}{\partial }^{\nu }\left(\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\mu \nu }^{\left(j\right)bc}(x,x)\right]+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\mu }^{\left(j\right)b}\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\nu }^{\left(j\right)c}\left(x\right)\right]\right)\phantom{(}\hfill \\ \phantom{(}+e^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}({\partial }_{\mu }{\underline{G}}_{2\nu b}^{\left(j\right)c\nu }(x,x)-{\partial }_{\nu }{\underline{G}}_{2\mu b}^{\left(j\right)c\nu }(x,x))\right]\phantom{(}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1b}^{\left(j\right)\nu }\left(x\right)\right]e^{-f(\square )}\left[e^{\frac{1}{2}f(\square )}({\partial }_{\mu }{\underline{G}}_{1\nu }^{\left(j\right)c}\left(x\right)-{\partial }_{\nu }{\underline{G}}_{1\mu }^{\left(j\right)c}\left(x\right))\right]\}\phantom{(}\hfill \\ \phantom{(}+g_s^2{f}_{abc}{f}_{cde}\left\{e^{-\frac{1}{2}f(\square )}\right(\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{3\mu \nu b}^{\left(j\right)de\nu }(x,x,x)\right]+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\mu b}^{\left(j\right)d\nu }(x,x)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\nu }^{\left(j\right)e}\left(x\right)\right]\phantom{(}\hfill \\ \phantom{(}+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\mu }^{\left(j\right)d}\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\nu b}^{\left(j\right)e\nu }(x,x)\right])\phantom{(}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1b}^{\left(j\right)\nu }\left(x\right)\right]e^{-f(\square )}(\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\mu \nu }^{\left(j\right)de}(x,x)\right]\phantom{(}\hfill \\ \phantom{(}+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\mu }^{\left(j\right)d}\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\nu }^{\left(j\right)e}\left(x\right)\right]\left)\right\}.\hfill \end{array}$$

(13)

The equation of motion for the two-point function is obtained by variation with respect to $j_h^{\lambda }\left(y\right)$:

$$\begin{array}{c}{\displaystyle \square {\underline{G}}_{2\mu \lambda }^{\left(j\right)ah}(x,y)-i{\delta }^{ah}{\eta }_{\mu \lambda }{e}^{\frac{1}{2}f(\square )}{\delta }^{\left(4\right)}(x-y)=}\hfill \\ =g_s{f}_{abc}\left\{e^{-\frac{1}{2}f(\square )}{\partial }^{\nu }\right(\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{3\mu \nu \lambda }^{\left(j\right)bch}(x,x,y)\right]+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\mu \lambda }^{\left(j\right)bh}(x,y)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\nu }^{\left(j\right)c}\left(x\right)\right]\phantom{(}\hfill \\ \phantom{(}+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\mu }^{\left(j\right)b}\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\nu \lambda }^{\left(j\right)ch}(x,y)\right])\phantom{(}\hfill \\ \phantom{(}+e^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}({\partial }_{\mu }{\underline{G}}_{3\nu b\lambda }^{\left(j\right)c\nu h}(x,x,y)-{\partial }_{\nu }{\underline{G}}_{3\mu b\lambda }^{\left(j\right)c\nu h}(x,x,y))\right]\phantom{(}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}(\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2b\lambda }^{\left(j\right)\nu h}(x,y)\right]e^{-f(\square )}\left[e^{\frac{1}{2}f(\square )}({\partial }_{\mu }{\underline{G}}_{1\nu }^{\left(j\right)c}\left(x\right)-{\partial }_{\nu }{\underline{G}}_{1\mu }^{\left(j\right)c}\left(x\right))\right]\phantom{(}\hfill \\ \phantom{(}+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1b}^{\left(j\right)\nu }\left(x\right)\right]e^{-f(\square )}\left[e^{\frac{1}{2}f(\square )}({\partial }_{\mu }{\underline{G}}_{2\nu \lambda }^{\left(j\right)ch}(x,y)-{\partial }_{\nu }{\underline{G}}_{2\mu \lambda }^{\left(j\right)ch}(x,y))\right]\left)\right\}\phantom{(}\hfill \\ \phantom{(}+g_s^2{f}_{abc}{f}_{cde}\left\{e^{-\frac{1}{2}f(\square )}\right(\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{4\mu \nu b\lambda }^{\left(j\right)de\nu h}(x,x,x,y)\right]\phantom{(}\hfill \\ \phantom{(}+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{3\mu \nu b}^{\left(j\right)d\nu h}(x,x,y)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\nu }^{\left(j\right)e}\left(x\right)\right]+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\mu b}^{\left(j\right)d\nu }(x,x)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\nu \lambda }^{\left(j\right)eh}(x,y)\right]\phantom{(}\hfill \\ \phantom{(}+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\mu \lambda }^{\left(j\right)dh}(x,y)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\nu b}^{\left(j\right)e\nu }(x,x)\right]+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\mu }^{\left(j\right)d}\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{3\nu b\lambda }^{\left(j\right)e\nu h}(x,x,y)\right])\phantom{(}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2b\lambda }^{\left(j\right)\nu h}(x,y)\right]\times \phantom{(}\hfill \\ \phantom{(}\times e^{-f(\square )}(\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\mu \nu }^{\left(j\right)de}(x,x)\right]+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\mu }^{\left(j\right)d}\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\nu }^{\left(j\right)e}\left(x\right)\right])\phantom{(}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1b}^{\left(j\right)\nu }\left(x\right)\right]e^{-f(\square )}(\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{3\mu \nu \lambda }^{\left(j\right)deh}(x,x,y)\right]\phantom{(}\hfill \\ \phantom{(}+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\mu \lambda }^{\left(j\right)dh}(x,y)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\nu }^{\left(j\right)e}\left(x\right)\right]+\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{1\mu }^{\left(j\right)d}\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_{2\nu \lambda }^{\left(j\right)eh}(x,y)\right]\left)\right\},\hfill \end{array}$$

(14)

with $\delta {\underline{G}}_{1\mu }^{\left(j\right)a}\left(x\right)/\delta j_h^{\lambda }\left(y\right)={\underline{G}}_{2\mu \lambda }^{\left(j\right)ah}(x,y)$. The system of equations of motion for the complete set of components of the Green functions for the Yang–Mills field cannot be treated exactly. Instead, we use a mapping to the scalar case. The mapping theorem introduced in Refs. [86,87] is based on Andrei Smilga’s solution for the problem that “the Yang–Mills system is not exactly solvable, … in contrast to … some early hopes” (cf. Sec. 1.2 in Ref. [88]). Indeed, for more than one independent component, in solving the system of equations, one observes chaotic behaviour (cf. Refs. [89,90,91]). Even though the importance of such observations for macroscopic behaviour is questionable, the safe path is to use a mapping to a single scalar function that we dub $\underline{\varphi }\left(x\right)$, representing an exact solution for the 1P-correlation function in the scalar case. This is realized here by using ${\underline{G}}_{1\mu }^{\left(0\right)a}\left(x\right)={\eta }_{\mu }^a\underline{\varphi }\left(x\right)$ and ${\underline{G}}_{2\mu \nu }^{\left(0\right)ab}(x,y)={\eta }_{\mu \nu }{\delta }^{ab}{\underline{G}}_2(x-y)$ for $j_{\mu }^a\left(x\right)$ set to zero, where ${\eta }_{\mu }^a$ represents the components of the polarisation vector and ${\eta }_{\mu \nu }$ denotes the components of the Minkowski metric. The two-point Green functions from x to x become constants with a vanishing derivative, while n-point Green functions with $n>2$ can be set to zero if at least two arguments coincide. One obtains

$$\begin{array}{c}{\displaystyle {\eta }_{\mu }^a\square \underline{\varphi }\left(x\right)=g_s{f}_{abc}\{e^{-\frac{1}{2}f(\square )}{\partial }^{\nu }{\eta }_{\mu }^b{\eta }_{\nu }^c{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\phantom{(}}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}{\eta }_b^{\nu }\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]e^{-f(\square )}\left[e^{\frac{1}{2}f(\square )}({\eta }_{\nu }^c{\partial }_{\mu }\underline{\varphi }\left(x\right)-{\eta }_{\mu }^c{\partial }_{\nu }\underline{\varphi }\left(x\right))\right]\}\phantom{(}\hfill \\ \phantom{(}+g_s^2{f}_{abc}{f}_{cde}\{e^{-\frac{1}{2}f(\square )}({\eta }_{\mu }^{\nu }{\delta }_b^d{\eta }_{\nu }^e+{\eta }_{\nu }^{\nu }{\delta }_b^e{\eta }_{\mu }^d)\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-x)\right]\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]\phantom{(}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}{\eta }_b^{\nu }\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]e^{-f(\square )}\hfill \\ \phantom{(}\times \left({\eta }_{\mu \nu }{\delta }^{de}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-x)\right]+{\eta }_{\mu }^d{\eta }_{\nu }^e{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\right)\}\hfill \\ =g_s{f}_{abc}\{e^{-\frac{1}{2}f(\square )}{\eta }_{\mu }^b{\eta }_c^{\nu }{\partial }_{\nu }{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\phantom{(}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]e^{-f(\square )}\left[e^{\frac{1}{2}f(\square )}({\eta }_b^{\nu }{\eta }_{\nu }^c{\partial }_{\mu }\underline{\varphi }\left(x\right)-{\eta }_b^{\nu }{\eta }_{\mu }^c{\partial }_{\nu }\underline{\varphi }\left(x\right))\right]\}\phantom{(}\hfill \\ \phantom{(}+g_s^2{f}_{abc}{f}_{cde}\{({\delta }^{bd}{\eta }_{\mu }^e+D{\delta }^{be}{\eta }_{\mu }^d)e^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-x)\right]\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]\phantom{(}\hfill \\ \phantom{(}+{\eta }_b^{\nu }{\eta }_{\mu }^d{\eta }_{\nu }^e{e}^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]e^{-f(\square )}{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\}\hfill \end{array}$$

(15)

and

$$\begin{array}{c}{\displaystyle {\delta }^{ah}{\eta }_{\mu \lambda }\square {\underline{G}}_2(x-y)-i{\delta }^{ah}{\eta }_{\mu \lambda }{e}^{\frac{1}{2}f(\square )}{\delta }^{\left(4\right)}(x-y)=}\hfill \\ =g_s{f}_{abc}\{e^{-\frac{1}{2}f(\square )}({\eta }_{\mu \lambda }{\delta }^{bh}{\eta }_{\nu }^c+{\eta }_{\nu \lambda }{\eta }_{\mu }^b{\delta }^{ch}){\partial }^{\nu }\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-y)\right]\phantom{(}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}{\delta }_b^h{\eta }_{\lambda }^{\nu }\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-y)\right]e^{-f(\square )}\left[e^{\frac{1}{2}f(\square )}({\eta }_{\nu }^c{\partial }_{\mu }\underline{\varphi }\left(x\right)-{\eta }_{\mu }^c{\partial }_{\nu }\underline{\varphi }\left(x\right))\right]\phantom{(}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}{\eta }_b^{\nu }{\delta }^{ch}\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]e^{-f(\square )}\left[e^{\frac{1}{2}f(\square )}({\eta }_{\nu \lambda }{\partial }_{\mu }-{\eta }_{\mu \lambda }{\partial }_{\nu }){\underline{G}}_2(x-y)\right]\phantom{(}\hfill \\ \phantom{(}+g_s^2{f}_{abc}{f}_{cde}\{{\eta }_{\mu \lambda }({\delta }_b^d{\delta }^{eh}+{\delta }_b^e{\delta }^{dh}{\eta }_{\nu }^{\nu })\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-x)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-y)\right]\phantom{(}\hfill \\ \phantom{(}+{\eta }_{\lambda }^{\nu }{\delta }_b^h{e}^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-y)\right]\hfill \\ \phantom{(}\times e^{-f(\square )}\left({\eta }_{\mu \nu }{\delta }^{de}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-x)\right]+{\eta }_{\mu }^d{\eta }_{\nu }^e{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\right)\phantom{(}\hfill \\ \phantom{(}+{\eta }_b^{\nu }({\eta }_{\mu \lambda }{\delta }^{dh}{\eta }_{\nu }^e+{\eta }_{\nu \lambda }{\delta }^{eh}{\eta }_{\mu }^d)\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-y)\right]{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\}.\hfill \end{array}$$

(16)

Contracting with ${\eta }_a^{\mu }$ and ${\delta }_{ah}{\eta }^{\mu \lambda }$ and using the orthogonality and completeness relations

$${\eta }_{\mu }^a{\eta }_b^{\mu }=-{\delta }_b^a,{\eta }_{\mu }^a{\eta }_{\nu }^a=-(N_c^2-1){\eta }_{\mu \nu }/D,$$

(17)

where $\left({\eta }_{\mu \nu }\right)=diag(1,-1,-1,-1)$ ($D={\eta }_{\mu \lambda }{\eta }^{\mu \lambda }$ is the space–time dimension), one obtains

$$\begin{array}{c}-(N_c^2-1)\square \underline{\varphi }\left(x\right)=g_s^2{f}_{abc}{f}_{cde}\{-({\delta }_b^d{\delta }_a^e+D{\delta }_b^e{\delta }_a^d)e^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-x)\right]\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]\phantom{(}\hfill \\ \phantom{(}+{\delta }_b^e{\delta }_a^d{e}^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]e^{-f(\square )}{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\}\hfill \\ =N_c(N_c^2-1)g_s^2\{-(D-1)e^{-\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-x)\right]\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]\phantom{(}\hfill \\ \phantom{(}+e^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]e^{-f(\square )}{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\}\hfill \end{array}$$

(18)

and

$$\begin{array}{c}{\displaystyle D(N_c^2-1)\left(\square {\underline{G}}_2(x-y)-i{e}^{\frac{1}{2}f(\square )}{\delta }^{\left(4\right)}(x-y)\right)=}\hfill \\ =g_s^2{f}_{abc}{f}_{cde}\{D({\delta }_b^d{\delta }_a^e+D{\delta }_b^e{\delta }_a^d)\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-x)\right]\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-y)\right]\phantom{(}\hfill \\ \phantom{(}-(D{\delta }_a^d{\delta }_b^e+{\delta }_a^e{\delta }_b^d)\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-y)\right]{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\}\hfill \\ =(D-1)N_c(N_c^2-1)g_s^2\left\{D\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-x)\right]-{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\right\}\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-y)\right],\hfill \end{array}$$

(19)

where $f_{abc}{f}_{abd}=N_c{\delta }_{cd}$ for two and, therefore, $f_{abc}{f}_{abc}=N_c(N_c^2-1)$ for three summed indices have been used. The common factor cancels, and one obtains

$$(\square +\Delta m_G^2)\underline{\varphi }\left(x\right)+\lambda e^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]e^{-f(\square )}{\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2=0$$

(20)

and

$$(\square +\Delta m_G^2){\underline{G}}_2(x-y)+\frac{D-1}{D}\lambda {\left[e^{\frac{1}{2}f(\square )}\underline{\varphi }\left(x\right)\right]}^2\left[e^{\frac{1}{2}f(\square )}{\underline{G}}_2(x-y)\right]=i{e}^{\frac{1}{2}f(\square )}{\delta }^{\left(4\right)}(x-y),$$

(21)

where $\lambda =N_c{g}_s^2$ and $\Delta m_G^2=-(D-1)e^{\frac{1}{2}f\left({\square }_x\right)}\lambda G_2(x-x^{\prime }){|}_{x^{\prime }=x}$.

Looking at Equation (20), the solution for the corresponding local one-point function, given by ${\underline{\varphi }}_0\left(x\right)={\widehat{\varphi }}_0\left(kx\right)=\mu \mathrm{sn}(kx+\theta |\kappa )$ [86], obeys the dispersion relation $k^2=\Delta m_G^2+\lambda {\mu }^2/2$, with $\kappa =(\Delta m_G^2-k^2)/k^2$ and $\theta =(1+4N)K\left(\kappa \right)$. The newly introduced constant $\mu$ needs some explanation. This is an integration constant coming from the solution of the Dyson–Schwinger equation for the one-point correlation function, with the dimension of energy. It sets the scale for our non-perturbative solution. Therefore, it should be related to the constant ${\Lambda }_{\mathrm{QCD}}$ that emerges from dimensional transmutation in perturbative computations, setting the scale for the confined phase of Yang–Mills theory. The solution ${\widehat{\varphi }}_0\left(kx\right)$ can be expanded in a Fourier series,

$$\begin{array}{ccc}\hfill {\widehat{\varphi }}_0\left(kx\right)& =& \frac{2\eta }{\sqrt{\kappa }}\sqrt{\frac{-2k^2\kappa }{\lambda }}\sum _{n=0}^{\infty }\frac{q^{n+1/2}}{1-q^{2n+1}}sin\left(\frac{2n+1}{2}\eta (kx+\theta )\right)\hfill \\ & =& -2i\eta \sqrt{\frac{2k^2}{\lambda }}\sum _{n=0}^{\infty }\frac{{(-1)}^n{q}^{n+1/2}}{1-q^{2n+1}}cos\left(\frac{2n+1}{2}\eta kx\right)\hfill \\ & =& -i\mu \eta \sum _{n=0}^{\infty }\frac{{(-1)}^n{q}^{n+1/2}}{1-q^{2n+1}}(e^{(n+1/2)i\eta kx}+e^{-(n+1/2)i\eta kx})\hfill \\ & =& -i\mu \eta \sum _{n=-\infty }^{\infty }{b}_{2n+1}{e}^{(n+1/2)i\eta kx}=-i\mu \eta \sum _{m\mathrm{odd}}{b}_m{e}^{im\eta kx/2},\hfill \end{array}$$

(22)

with $\eta =\pi /K\left(\kappa \right)$. For non-negative values of n, one obviously has

$$b_{2n+1}=\frac{{(-1)}^n{q}^{n+1/2}}{1-q^{2n+1}}.$$

(23)

However, this can also be extended to negative values, such as for

$$\begin{array}{ccc}\hfill b_{-2n-1}& =& \frac{{(-1)}^{n-1}{q}^{-n-1/2}}{1-q^{-2n-1}}=\frac{-{(-1)}^n{q}^{2n+1}{q}^{-n-1/2}}{q^{2n+1}(1-q^{-2n-1})}\hfill \\ & =& \frac{-{(-1)}^n{q}^{n+1/2}}{q^{2n+1}-1}=\frac{{(-1)}^n{q}^{n+1/2}}{1-q^{2n+1}}=b_{2n+1}.\hfill \end{array}$$

(24)

First-order non-locality effects are taken into account by solving

$$(\square +\Delta m_G^2){\underline{\varphi }}_1\left(x\right)=-\lambda e^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{\varphi }}_0\left(x\right)\right]e^{-f(\square )}{\left[e^{\frac{1}{2}f(\square )}{\underline{\varphi }}_0\left(x\right)\right]}^2.$$

(25)

For the solution ${\underline{\varphi }}_1\left(x\right)={\widehat{\varphi }}_1\left(kx\right)$, we can use a similar ansatz:

$${\widehat{\varphi }}_1\left(kx\right)=-i\eta \mu \sum _{m\mathrm{odd}}{b}_m{e}^{im\eta kx/2},$$

(26)

where the $b_m$ values are to be computed. Substituting ${\underline{\varphi }}_0\left(x\right)={\widehat{\varphi }}_0\left(kx\right)$ on the right-hand side and ${\underline{\varphi }}_1\left(x\right)={\widehat{\varphi }}_1\left(kx\right)$ on the left-hand side, on the right-hand side, one calculates, step by step,

$$\begin{array}{c}{\displaystyle {\left[e^{\frac{1}{2}f(\square )}{\underline{\varphi }}_0\left(x\right)\right]}^2=-{\eta }^2{\mu }^2\sum _{m_2,m_3\mathrm{odd}}{b}_{m_2}{b}_{m_3}{e}^{\frac{1}{2}f(-m_2^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_3^2{\eta }^2{k}^2/4)}{e}^{i(m_2+m_3)\eta kx/2},}\hfill \\ {\displaystyle e^{-f(\square )}{\left[e^{\frac{1}{2}f(\square )}{\underline{\varphi }}_0\left(x\right)\right]}^2=-{\eta }^2{\mu }^2\sum _{m_2,m_3\mathrm{odd}}{b}_{m_2}{b}_{m_3}{e}^{-f(-{(m_2+m_3)}^2{\eta }^2{k}^2/4)}\phantom{(}}\hfill \\ \phantom{(}\times e^{\frac{1}{2}f(-m_2^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_3^2{\eta }^2{k}^2/4)}{e}^{i(m_2+m_3)\eta kx/2},\hfill \\ {\displaystyle \left[e^{\frac{1}{2}f(\square )}{\underline{\varphi }}_0\left(x\right)\right]e^{-f(\square )}{\left[e^{\frac{1}{2}f(\square )}{\underline{\varphi }}_0\left(x\right)\right]}^2=i{\eta }^3{\mu }^3\sum _{m_1,m_2,m_3\mathrm{odd}}{b}_{m_1}{b}_{m_2}{b}_{m_3}\phantom{(}}\hfill \\ \phantom{(}\times e^{\frac{1}{2}f(-m_1^2{\eta }^2{k}^2/4)}{e}^{-f(-{(m_2+m_3)}^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_2^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_3^2{\eta }^2{k}^2/4)}{e}^{i(m_1+m_2+m_3)\eta kx/2},\hfill \\ {\displaystyle e^{\frac{1}{2}f(\square )}\left[e^{\frac{1}{2}f(\square )}{\underline{\varphi }}_0\left(x\right)\right]e^{-f(\square )}{\left[e^{\frac{1}{2}f(\square )}{\underline{\varphi }}_0\left(x\right)\right]}^2}\hfill \\ =i{\eta }^3{\mu }^3\sum _{m_1,m_2,m_3\mathrm{odd}}{b}_{m_1}{b}_{m_2}{b}_{m_3}{e}^{\frac{1}{2}f(-{(m_1+m_2+m_3)}^2{\eta }^2{k}^2/4)}\phantom{(}\hfill \\ \phantom{(}\times e^{\frac{1}{2}f(-m_1^2{\eta }^2{k}^2/4)}{e}^{-f(-{(m_2+m_3)}^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_2^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_3^2{\eta }^2{k}^2/4)}{e}^{i(m_1+m_2+m_3)\eta kx/2}.\hfill \end{array}$$

(27)

Therefore, the coefficients $b_m$ are determined by the following system of equations:

$$\begin{array}{c}{\displaystyle b_m(-m^2{\eta }^2{k}^2/4+\Delta m_G^2)=\lambda {\eta }^2{\mu }^2\sum _{m_1,m_2,m_3\mathrm{odd}}^{m_1+m_2+m_3=m}{b}_{m_1}{b}_{m_2}{b}_{m_3}\phantom{(}}\hfill \\ \phantom{(}\times e^{\frac{1}{2}f(-{(m_1+m_2+m_3)}^2{\eta }^2{k}^2/4)}{e}^{-f(-{(m_2+m_3)}^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_1^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_2^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_3^2{\eta }^2{k}^2/4)}.\hfill \end{array}$$

(28)

One can solve for $b_m$ to obtain

$$\begin{array}{c}{\displaystyle b_m=\lambda {\eta }^2{\mu }^2\sum _{m_1,m_2,m_3\mathrm{odd}}^{m_1+m_2+m_3=m}{b}_{m_1}{b}_{m_2}{b}_{m_3}{(-{(m_1+m_2+m_3)}^2{\eta }^2{k}^2/4+\Delta m_G^2)}^{-1}\phantom{(}}\hfill \\ \phantom{(}\times e^{\frac{1}{2}f(-{(m_1+m_2+m_3)}^2{\eta }^2{k}^2/4)}{e}^{-f(-{(m_2+m_3)}^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_1^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_2^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_3^2{\eta }^2{k}^2/4)},\hfill \end{array}$$

(29)

and after inserting it back into ${\underline{\varphi }}_1\left(x\right)$, as the first iteration, one has

$$\begin{array}{c}{\displaystyle {\underline{\varphi }}_1\left(x\right)=i\lambda {\eta }^3{\mu }^3\sum _{m_1,m_2,m_3\mathrm{odd}}\frac{b_{m_1}{b}_{m_2}{b}_{m_3}{e}^{i(m_1+m_2+m_3)\eta kx/2}}{{(m_1+m_2+m_3)}^2{\eta }^2{k}^2/4-\Delta m_G^2}\phantom{(}}\hfill \\ \phantom{(}\times e^{\frac{1}{2}f(-{(m_1+m_2+m_3)}^2{\eta }^2{k}^2/4)}{e}^{-f(-{(m_2+m_3)}^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_1^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_2^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_3^2{\eta }^2{k}^2/4)}.\hfill \end{array}$$

(30)

However, this is not the physical field. Instead, going back to physics, we have to revert the field redefinition, obtaining ${\varphi }_1\left(x\right)=e^{\frac{1}{2}f(\square )}{\underline{\varphi }}_1\left(x\right)$ or

$$\begin{array}{c}{\displaystyle {\varphi }_1\left(x\right)=i\lambda {\eta }^3{\mu }^3\sum _{m_1,m_2,m_3\mathrm{odd}}\frac{b_{m_1}{b}_{m_2}{b}_{m_3}{e}^{i(m_1+m_2+m_3)\eta kx/2}}{{(m_1+m_2+m_3)}^2{\eta }^2{k}^2/4-\Delta m_G^2}\phantom{(}}\hfill \\ \phantom{(}\times e^{f(-{(m_1+m_2+m_3)}^2{\eta }^2{k}^2/4)}{e}^{-f(-{(m_2+m_3)}^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_1^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_2^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-m_3^2{\eta }^2{k}^2/4)}.\hfill \end{array}$$

(31)

Due to the denominator factor, for $f(\square )=\square /M^2$, the series turns out to be convergent.

The same is performed for Equation (21) in the case of $D=4$. After inserting $\underline{\varphi }\left(x\right)={\underline{\varphi }}_0\left(x\right)$, in momentum space, one has

$$\begin{array}{c}{\displaystyle -(p^2-\Delta m_G^2{e}^{\frac{1}{2}f(-p^2)}){\underline{\tilde{G}}}_2\left(p\right)-i{e}^{\frac{1}{2}f(-p^2)}}\hfill \\ \phantom{(}=\frac{3}{4}\lambda {\mu }^2{\eta }^2\sum _{m,n\mathrm{odd}}{b}_m{b}_n{e}^{\frac{1}{2}f(-m^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-n^2{\eta }^2{k}^2/4)}{e}^{\frac{1}{2}f(-{(p+\frac{1}{2}(m+n)\eta k)}^2)}{\underline{\tilde{G}}}_2(p+\frac{1}{2}(m+n)\eta k).\hfill \end{array}$$

(32)

The localised equation is solved in Appendix B and, according to the Källén–Lehmann theorem, provides a spectrum of harmonics, known as glueball states, which can be understood as a form of the glueball operators in Yang–Mills theory according to our solution. Gaussian factors can be used to restrict the solution to the lowest harmonics $m,n=\pm 1$, resulting in

$$\begin{array}{c}{\displaystyle -(p^2-\Delta m_G^2{e}^{\frac{1}{2}f(-p^2)}){\underline{\tilde{G}}}_2\left(p\right)-i{e}^{\frac{1}{2}f(-p^2)}}\hfill \\ =-m_G^2\left(e^{\frac{1}{2}f(-p^2)}{\underline{\tilde{G}}}_2\left(p\right)+\frac{1}{2}{e}^{\frac{1}{2}f(-{(p-\eta k)}^2)}{\underline{\tilde{G}}}_2(p-\eta k)+\frac{1}{2}{e}^{\frac{1}{2}f(-{(p+\eta k)}^2)}{\underline{\tilde{G}}}_2(p+\eta k)\right),\hfill \end{array}$$

(33)

where

$$m_G^2=-\frac{3}{2}\lambda {\mu }^2{\eta }^2{b}_1^2{e}^{f(-{\eta }^2{k}^2/4)}=3\kappa {\eta }^2{k}^2{b}_1^2{e}^{f(-{\eta }^2{k}^2/4)}$$

(34)

The equation can be solved for ${\underline{\tilde{G}}}_2\left(p\right)$, resulting in

$${\underline{\tilde{G}}}_2\left(p\right)=\underline{\tilde{\Delta }}\left(p\right)\left[1+\frac{i}{2}{m}_G^2{e}^{-\frac{1}{2}f(-p^2)}\left(e^{\frac{1}{2}f(-{(p-\eta k)}^2)}{\underline{\tilde{G}}}_2(p-\eta k)+e^{\frac{1}{2}f(-{(p+\eta k)}^2)}{\underline{\tilde{G}}}_2(p+\eta k)\right)\right],$$

(35)

where

$$\underline{\tilde{\Delta }}\left(p\right):=\frac{-i{e}^{\frac{1}{2}f(-p^2)}}{p^2-(m_G^2+\Delta m_G^2)e^{\frac{1}{2}f(-p^2)}}.$$

(36)

Equation (35) can be solved iteratively up to arbitrary orders of $m_G^2$, as shown in Ref. [56]. For simplicity, we take only the leading-order approximation ${\underline{\tilde{G}}}_2\left(p\right)\approx \underline{\tilde{\Delta }}\left(p\right)$.

4. The Mass Gap Equation for the Glueball

The mass gap equation obtained from $\Delta m_G^2=-(D-1)\lambda e^{\frac{1}{2}f\left({\square }_x\right)}{G}_2(x-x^{\prime }){|}_{x^{\prime }=x}$ reads

$$\begin{array}{ccc}\hfill \Delta m_G^2& =& \int \frac{d^{4}p}{{\left(2\pi \right)}^4}\frac{3i\lambda e^{f(-p^2)}}{p^2-(m_G^2+\Delta m_G^2)e^{\frac{1}{2}f(-p^2)}}=-i\int \frac{d^4{p}_E}{{\left(2\pi \right)}^4}\frac{3i\lambda e^{f\left(p_E^2\right)}}{-p_E^2-(m_G^2+∆m_G^2)e^{\frac{1}{2}f\left(p_E^2\right)}}\hfill \\ & =& \frac{-3\lambda }{{\left(4\pi \right)}^2}{\int }_0^{\infty }\frac{{\rho }^3{e}^{f\left({\rho }^2\right)}d\rho }{{\rho }^2+(m_G^2+\Delta m_G^2)e^{\frac{1}{2}f\left({\rho }^2\right)}}.\hfill \end{array}$$

(37)

Note that both $m_G^2$ as defined before and $\Delta m_G^2=(1+\kappa )k^2$ contain the common factor $k^2$. This common factor sets the scale of our solution, as it is proportional to ${\mu }^2$, the scale of confinement in Yang–Mills theory, as discussed above. Therefore, if we were to also consider quarks, $\mu$ would be related to ${\Lambda }_{\mathrm{QCD}}$ as the scale of asymptotic freedom. The lattice community will recognise this common factor as what they call the “string tension”. At this stage, we do not want to evaluate the inter-quark potential to see the full relation between all three of these experimental constants. What we need to compare with lattice data are just the pure numbers arising from the ratio with such general confinement scales. While the integral is regular at the lower limit, it has to be regularised for the upper limit. We emphasise that the fact that there is a non-local factor does not imply at all that all the integrals in the theory are UV-finite. What such an approach grants is a physical cut-off scale, M, given by the theory, working like a horizon to keep such integrals finite. Dimensional regularisation is not applicable in this case, as the integral cannot be computed analytically, in contrast to the local limit $M\to \infty$. One might think of ${\rho }^2=M^2$ as an upper cut-off. In this case, the integral will diverge as $M^2$, and a truncation of this highest power is necessary in order to obtain a finite local limit. Instead, our proposal is to use an upper cut-off ${\Lambda }^2$ fixed up to the scale $k^2$, which is determined by matching the results to values obtained on the lattice for different values of $N_c=2,3,4,5,6,8,10,12$ [92]. Results from an older Ref. [93] were analysed in the local case in Refs. [94,95], with excellent agreement. The idea is to fix the cut-off scale for the mass gap integral using the results of the local theory (given for $M\to \infty$) in order to prevent changes to the physics with the introduction of the non-locality. Such changes, if any, do not appear experimentally, and they should be really tiny. In order to compare with the lattice, instead of $\lambda$ as the input, one has to use $\lambda =N_c{g}_s^2=2N_c^2/\beta$ with $\beta =2N_c/g_s^2$. The values for $\beta$ are taken from Ref. [92] and listed in the second column of

Table 1. Values for the dynamical glueball masses $M_G^{\mathrm{lat}}$ on the lattice in units of $\sqrt{k^2}$ for different values of the number $N_c$ of colours (and corresponding values of $\beta$), compared to our estimates $M_G^{\mathrm{est}}$ from the solution of the mass gap Equation (37) in the local limit $M\to \infty$.

${\mathit{N}}_{\mathit{c}}$	$\mathit{\beta }$	${\mathit{M}}_{\mathit{G}}^{\mathbf{lat}}/\sqrt{{\mathit{k}}^2}$	${\mathit{M}}_{\mathit{G}}^{\mathbf{est}}/\sqrt{{\mathit{k}}^2}$	Error
2	$2.427$	$3.781\left(23\right)$	$3.56137$	$-5.8\%$
3	$6.235$	$3.405\left(21\right)$	$3.25530$	$-4.4\%$
4	$11.02$	$3.271\left(27\right)$	$3.26827$	$-0.1\%$
5	$17.61$	$3.156\left(31\right)$	$3.21857$	$+2.0\%$
6	$25.35$	$3.102\left(32\right)$	$3.21930$	$+3.8\%$
8	$45.50$	$3.099\left(26\right)$	$3.19833$	$+3.2\%$
10	$71.38$	$3.102\left(37\right)$	$3.18957$	$+2.8\%$
12	$103.03$	$3.156\left(33\right)$	$3.18443$	$+1.1\%$

.

(A recent paper about the Casimir effect in non-Abelian gauge theories on the lattice [96] has shown that the ground state of local Yang–Mills theory is not the same as that found in lattice computations (cf. Table 1). The authors of Ref. [96] obtained $M_G^{\mathrm{lat}}/\sqrt{k^2}=1.0\left(1\right)$, which, for us, corresponds to the choice $n=0$ in the glueball spectrum displayed in Appendix B and Refs. [94,95]. This is also in agreement with the picture of the Casimir effect in local Yang–Mills theory presented in Ref. [97]). The agreement is excellent. For the fit of the mass values in the local limit, we used the method of least squares by minimising

$$\sigma \left(\Lambda \right)=\frac{1}{N}\sum _{i=1}^N{(M_G^{\left(i\right)\mathrm{est}}\left(\Lambda \right)-M_G^{\left(i\right)\mathrm{lat}})}^2,$$

(38)

with $M_G^2=m_G^2+\Delta m_G^2$. The result is given by ${\Lambda }^2=177.16\left(3\right)k^2$. As this result suggests, the agreement between the lattice values and our estimates is very good, as these are found in the remaining lines of Table 1. In Figure 1, we show the solution of the mass gap Equation (37) for $M_G^2$ for the upper limit ${\Lambda }^2=177.16k^2$ in relation to the non-local scale $M^2$, together with the dependence of the parameter $\kappa$. This could possibly be of the order of ${\mathrm{TeV}}^2$. It is interesting to point out that, in this way, we can obtain a physical understanding of the non-local scale and its proper order of magnitude as compared to local physics.

5. The Mass gap Equation of the Quark

The solution obtained in the previous section is input for the determination of the dynamical quark mass. Namely, the Green function ${\underline{G}}_{2\mu \nu }^{\left(0\right)ab}(x,y)={\eta }_{\mu \nu }{\delta }^{ab}{\underline{G}}_2(x-y)$, with ${\underline{G}}_2(x-y)$ satisfying the differential Equation (21), has to be convoluted with the left-hand side of the equation of motion (10) in order to obtain

$${\underline{A}}_a^{\nu }\left(x\right)=i{g}_s\int {\underline{G}}_2(x-y)e^{\frac{1}{2}f\left({\square }_y\right)}\overline{\psi }\left(y\right){\gamma }^{\nu }{T}_a\psi \left(y\right)d^{4}y.$$

(39)

This result is inserted into the Euler–Lagrange equation for the quark field,

$$0=\frac{\partial \mathcal{L}}{\partial \overline{\psi }}-{\partial }_{\mu }\frac{\partial \mathcal{L}}{\partial \left({\partial }_{\mu }\overline{\psi }\right)}=\frac{\partial \mathcal{L}}{\partial \overline{\psi }}={\gamma }^{\mu }\left(i{\partial }_{\mu }-g_s{T}_a{e}^{\frac{1}{2}f(\square )}{\underline{A}}_{\mu }^a\right)\psi$$

(40)

to obtain

$$0=i{\gamma }^{\mu }{\partial }_{\mu }\psi \left(x\right)+i{g}_s^2{\gamma }^{\mu }{T}_a\psi \left(x\right)e^{\frac{1}{2}f(\square )}\int {\underline{G}}_2(x-y)e^{\frac{1}{2}f\left({\square }_y\right)}\overline{\psi }\left(y\right){\gamma }_{\mu }{T}_a\psi \left(y\right)d^{4}y.$$

(41)

This equation of motion can be understood as the equation of motion of a non-local Nambu–Jona-Lasinio (NJL) model. After integration by parts, we derive the NJL Lagrangian

$${\mathcal{L}}_{\mathrm{NJL}}=\overline{\psi }\left(x\right)i{\gamma }^{\mu }{\partial }_{\mu }\psi \left(x\right)+i{g}_s^2\overline{\psi }\left(x\right){\gamma }^{\mu }{T}_a\psi \left(x\right)e^{\frac{1}{2}f(\square )}\int \left[e^{\frac{1}{2}f\left({\square }_y\right)}{\underline{G}}_2(x-y)\right]\overline{\psi }\left(y\right){\gamma }_{\mu }{T}_a\psi \left(y\right)d^{4}y.$$

(42)

The standard procedure to arrive at the mass gap equation is taken from Ref. [98] and consists of several steps. First of all, a Fierz rearrangement leads to the action integral $S[\psi ,\overline{\psi }]=S_0[\psi ,\overline{\psi }]+S_{\mathrm{int}}[\psi ,\overline{\psi }]$ with

$$\begin{array}{ccc}\hfill S_0[\psi ,\overline{\psi }]& =& \int d^{4}x{\overline{\psi }}^i\left(x\right)i{\gamma }^{\mu }{\partial }_{\mu }{\psi }^i\left(x\right)and\hfill \\ \hfill S_{\mathrm{int}}[\psi ,\overline{\psi }]& =& i{g}_s^2\int {\underline{G}}_2^f(x-y){\overline{\psi }}^i\left(x\right){\Gamma }_{\alpha }^{ij}{\psi }^j\left(y\right){\overline{\psi }}^j\left(y\right){\Gamma }_{ji}^{\alpha }{\psi }^i\left(x\right)d^{4}y\hfill \\ & =& i{g}_s^2\int d^{4}w\int d^{4}z\overline{\psi }(w+\frac{z}{2}){\Gamma }_{\alpha }\psi (w-\frac{z}{2}){\underline{G}}_2^f\left(z\right)\overline{\psi }(w-\frac{z}{2}){\Gamma }^{\alpha }\psi (w+\frac{z}{2})\hfill \end{array}$$

(43)

with ${\underline{G}}_2^f\left(z\right)=e^{f(\square )}{\underline{G}}_2\left(z\right)$. This action integral is an ingredient for the functional integral

$$\mathcal{Z}=\int \mathcal{D}\psi \mathcal{D}\overline{\psi }exp(-iS[\psi ,\overline{\psi }]).$$

(44)

In the following, we use $i{\underline{G}}_2^f\left(z\right)=G{\underline{C}}_2\left(z\right)/2$ with $\int {\underline{C}}_2\left(z\right)d^{4}z=1$. The quartic interaction term can be removed by introducing the meson field $\left({\varphi }^{\alpha }\right)=(\sigma ,\overrightarrow{\pi })$ via the factor

$$\int \mathcal{D}\varphi exp\left(-\frac{i}{2G}\int d^{4}z{\underline{C}}_2\left(z\right)\int d^{4}w{\varphi }_{\alpha }^{*}\left(w\right){\varphi }^{\alpha }\left(w\right)\right),$$

(45)

into the functional integral, interchanging the integrations over w and z and performing the functional “shift”:

$${\varphi }_{\alpha }\left(w\right)\to {\varphi }_{\alpha }\left(w\right)+g_{s}G\overline{\psi }(w-\frac{z}{2}){\Gamma }_{\alpha }\psi (w+\frac{z}{2}).$$

(46)

In doing so and returning in part to x and y, one ends up with the action functional

$$\begin{array}{c}{\displaystyle S[\psi ,\overline{\psi },\varphi ]=\frac{1}{2G}\int d^{4}z{\underline{C}}_2\left(z\right)\int d^{4}w{\varphi }_{\alpha }^{*}\left(w\right){\varphi }^{\alpha }\left(w\right)\phantom{(}}\hfill \\ \phantom{(}+\int d^{4}x\int d^{4}y\overline{\psi }\left(x\right)\left[{\delta }^{\left(4\right)}(x-y)i{\gamma }^{\mu }{\partial }_{\mu }-g_s{\underline{C}}_2(x-y)\mathrm{Re}{\varphi }_{\alpha }\left(\frac{x+y}{2}\right){\Gamma }^{\alpha }\right]\psi \left(y\right),\hfill \end{array}$$

(47)

The functional integral reads

$$\mathcal{Z}=\int \mathcal{D}\psi \mathcal{D}\overline{\psi }\mathcal{D}\varphi exp(-iS[\psi ,\overline{\psi },\varphi ]).$$

(48)

Performing a Fourier transform and integrating out the fermionic degrees of freedom, one ends up with the functional determinant and the bosonised action

$$\begin{array}{c}{\displaystyle S_{\mathrm{bos}}\left[\varphi \right]=\frac{1}{2G}\int \frac{d^{4}q}{{\left(2\pi \right)}^4}{\varphi }_{\alpha }^{*}\left(q\right){\varphi }^{\alpha }\left(q\right)\phantom{(}}\hfill \\ +ilndet\left[{\left(2\pi \right)}^4{\delta }^{\left(4\right)}(p^{\prime }-p){\gamma }^{\mu }{p}_{\mu }-\frac{1}{2}{g}_s{\underline{\tilde{C}}}_2\left(\frac{p^{\prime }+p}{2}\right)({\varphi }_{\alpha }(p^{\prime }-p)+{\varphi }_{\alpha }^{*}(p-p^{\prime })){\Gamma }^{\alpha }\right],\hfill \end{array}$$

(49)

where the logarithm of the determinant is understood to be between states $\langle p^{\prime }|\cdots |p\rangle$ in momentum space. While in the mean field approximation, the first term gives ${\overline{\sigma }}^2/2G$, the determinant (for $N_f=2$ flavours) leads to

$$det({\gamma }^{\mu }{p}_{\mu }-M_q\left(p\right)1\mathrm{l})={(p^2-M_q^2\left(p\right))}^{4N_c},$$

(50)

with $M_q\left(p\right)=g_s{\underline{\tilde{C}}}_2\left(p\right)\overline{\sigma }$. The variation in the action integral with respect to the mean value $\overline{\sigma }$ leads to

$$\overline{\sigma }=-8i{N}_c{g}_{s}G\int \frac{d^{4}p}{{\left(2\pi \right)}^4}\frac{{\underline{\tilde{C}}}_2\left(p\right)M_q\left(p\right)}{p^2-M_q^2\left(p\right)},$$

(51)

and the insertion into the definition of $M_q\left(p\right)$ finally leads to the mass gap equation

$$\begin{array}{ccc}\hfill M_q\left(p\right)& =& -8i{N}_c{g}_s^{2}G{\underline{\tilde{C}}}_2\left(p\right)\int \frac{d^{4}q}{{\left(2\pi \right)}^4}\frac{{\underline{\tilde{C}}}_2\left(q\right)M_q\left(q\right)}{q^2-M_q^2\left(q\right)}\hfill \\ & =& 16\lambda \frac{{\tilde{G}}_2\left(p\right)}{{\underline{\tilde{G}}}_2\left(0\right)}\int \frac{d^{4}q}{{\left(2\pi \right)}^4}\frac{{\underline{\tilde{G}}}_2\left(q\right)M_q\left(q\right)}{q^2-M_q^2\left(q\right)},\hfill \end{array}$$

(52)

where we insert ${\underline{\tilde{C}}}_2\left(p\right)=-2i{\underline{\tilde{G}}}_2\left(p\right)/G$ and use

$${\underline{G}}_2^f\left(z\right)=\int \frac{d^{4}q}{{\left(2\pi \right)}^4}{e}^{f(\square )}{\underline{\tilde{G}}}_2\left(q\right)e^{-iqz}$$

(53)

to obtain

$$\frac{G}{2}=i\int {\underline{G}}_2^f\left(z\right)d^{4}z=i\int \int \frac{d^{4}q}{{\left(2\pi \right)}^4}{e}^{f(-q^2)}{\underline{\tilde{G}}}_2\left(q\right)e^{-iqz}{d}^{4}z=i{\underline{\tilde{G}}}_2\left(0\right).$$

(54)

In order to solve the gap equation for the quark, we assume that the mass gap $M_q\left(p\right)$ does not depend explicitly on the momentum. In this case, one has to solve the equation

$$M_q=16\lambda \int \frac{d^{4}q}{{\left(2\pi \right)}^4}\frac{e^{f(-q^2)}{\underline{\tilde{G}}}_2\left(q\right)M_q}{q^2-M_q^2}.$$

(55)

In principle, $M_q$ can be cancelled on both sides, leading to a gap equation similar to Equation (3.7) in the original NJL publication [99]. However, in order to solve this equation iteratively, it is more appropriate to instead multiply it with $M_q$ in order to determine $M_q^2$ from

$$\begin{array}{ccc}\hfill M_q^2& =& 16N_c{g}_s^2\int \frac{d^{4}q}{{\left(2\pi \right)}^4}\frac{e^{f(-q^2)}{\underline{\tilde{G}}}_2\left(q\right)M_q^2}{q^2-M_q^2}=\int \frac{d^{4}q}{{\left(2\pi \right)}^4}\frac{-16i\lambda M_q^2{e}^{\frac{3}{2}f(-q^2)}}{(q^2-M_G^2{e}^{\frac{1}{2}f(-q^2)})(q^2-M_q^2)}\hfill \\ & =& \int \frac{d^4{q}_E}{{\left(2\pi \right)}^4}\frac{16\lambda M_q^2{e}^{\frac{3}{2}f\left(q_E^2\right)}}{(-q_E^2-M_G^2{e}^{\frac{1}{2}f\left(q_E^2\right)})(-q_E^2-M_q^2)}\hfill \\ & =& \frac{8\lambda }{{\left(4\pi \right)}^2}{\int }_0^{{\Lambda }^2}\frac{M_q^2{\rho }^2{e}^{\frac{3}{2}f\left({\rho }^2\right)}d{\rho }^2}{({\rho }^2+M_G^2{e}^{\frac{1}{2}f\left({\rho }^2\right)})({\rho }^2+M_q^2)}.\hfill \end{array}$$

(56)

The upper limit is taken to be the same as in the previous section, and the solution $M_G^2$ of the mass gap equation for the glueball is used. The fixed-point problem converges, and one obtains a dependence on the non-local scale $M^2$, which is displayed in Figure 2. We can see that chiral symmetry breaking appears in non-local QCD as well.

It is shown in Refs. [58,59] that an educated guess for confinement is given by the threshold $M_q/M_G>0.39$, above which two real pole solutions for $p^2$ in the equation

$$p^2=M_q^2{\left(\frac{{\tilde{G}}_2\left(p\right)}{{\tilde{G}}_2\left(0\right)}\right)}^2=\frac{M_q^2{M}_G^4}{{(p^2-M_G^2)}^2}$$

(57)

change to imaginary pole solutions. According to Figure 3, this is the case for the non-local scale $M^2$ below $150k^2$.

6. Discussion and Conclusions

The spectrum of the Yang–Mills field is analysed for the non-local Yang–Mills sector and found to be in agreement with the local results on the lattice in the limit of the non-local scale running to infinity (cf. Table 1). The agreement with the lattice is astonishingly good, as is the agreement with previous computations [94,95]. Because of this, this spectrum becomes our reference for the non-local case. We point out that recent lattice results on the Casimir effect in Yang–Mills theory seem to challenge the determination of the ground state of the theory on the lattice. In our case, this would be a further confirmation of our approach (for details, see Ref. [96]). It is a relevant result of our investigations that scale invariance is badly broken by interactions, a fact that can be taken as a possible clue of confinement. Indeed, the solution of the gap equation for the fermion shows some indication of quark confinement in the non-local case as well. This result is really important, as it seems to point to the fact that confinement could be a ubiquitous effect in nature that removes degrees of freedom in a theory to favour others. A rigorous mathematical proof of confinement is beyond the scope of the present manuscript. In general, further studies are needed to improve these results. Still, the results obtained here appear to be a sound confirmation of previous work with a different technique and with the important extension to non-local QCD.