1 Introduction

The Faddeev–Popov operator \(\mathcal {M}\equiv \mathcal {M}^{ab}[A]=-\partial _\mu D_\mu ^{ab}\), with \(D_\mu ^{ab}=\delta ^{ab}\partial _\mu -gf^{abc}A_\mu ^c\) the covariant derivative in the adjoint representation, plays an important role when discussing geometric features of the Landau gauge fixing \(\partial _\mu A_\mu ^a=0\) for non-Abelian gauge theories. The set of all \(A_\mu \equiv A_\mu ^a\) that fulfill

$$\begin{aligned} \partial _\mu A_\mu =0 ~\wedge ~ \mathcal {M}\ge 0 \end{aligned}$$
(1)

defines the so-called Gribov region \(\Omega \). At its boundary \(\partial \Omega \), the Faddeev–Popov operator becomes zero. The definition (1) makes sense, since \(\mathcal {M}\) is a Hermitian operator, thus having real eigenvalues. It is then a simple matter to see that \(A_\mu ^a\in \Omega \) is a necessary and sufficient condition to ensure that \(A_\mu \) has no infinitesimally gauge equivalent field configurations also subject to the Landau gauge. The Gribov gauge fixing ambiguity at the infinitesimal level is thus overcome if the set of admissible gauge fields in the path integral are those belonging to \(\Omega \). A few interesting properties of the region \(\Omega \) have been established: it is convex, bounded in every direction (it is a kind of multidimensional ellipsoid) and most importantly, each gauge orbit intersects \(\Omega \) at least once [1, 2].

The conditions (1) are nothing else than those for identifying a local minimum of the functional [3, 4]

$$\begin{aligned} \mathcal {R}[A]=\min _{U\in \text {SU}(N)}\int \mathrm {d}^4 x A_\mu ^U A_\mu ^U. \end{aligned}$$
(2)

If we think in terms of the action of Landau gauge-fixed Yang–Mills theory,

$$\begin{aligned}&S_{\mathrm{GF}+\mathrm{FP}} \nonumber \\&\quad = \frac{1}{4}\int \mathrm {d}^d x\; F^a_{\mu \nu } F^a_{\mu \nu }+\int \mathrm {d}^d x\,\left( b^a \partial _\mu A_\mu ^a +\overline{c}^a \partial _\mu D_\mu ^{ab} c^b \right) ,\nonumber \\ \end{aligned}$$
(3)

the inverse Faddeev–Popov operator can be identified with the ghost propagator \(\mathcal {G}^{ab}(k,A)\) in the presence of an external gauge field. This observation is the crux of the original recipe [5] of Gribov to work out a practical version of the restriction to \(\Omega \): only gauge fields with \(\sigma (k,A)\le 1\) must be taken into account in the path integral, where \(G(k,A)=\frac{1}{N^2-1}\mathcal {G}^{aa}(k,A)=\frac{1}{k^2}\frac{1}{1-\sigma (k,A)}\). This is known as the no-pole condition and it was translated into an effective action at tree level by Gribov. Using an alternative derivation, Zwanziger [6, 7] was able to generalize the Gribov effective action to all orders. The end point was the following (Gribov–Zwanziger) action:

$$\begin{aligned} S_\mathrm {GZ}&= \frac{1}{4}\int \mathrm {d}^d x\; F^a_{\mu \nu } F^a_{\mu \nu }\!+\!\int \mathrm {d}^d x\,\left( b^a \partial _\mu A_\mu ^a\! +\!\overline{c}^a \partial _\mu D_\mu ^{ab} c^b \right) \nonumber \\&+\int \mathrm {d}^d x \left( \overline{\varphi }_{\mu }^{ac} \partial _{\nu } D_\nu ^{am}\varphi _{\mu }^{mc} -\overline{\omega }_{\mu }^{ac} \partial _{\nu } D_\nu ^{am} \omega _{\mu }^{mc} \right. \nonumber \\&-\left. g\left( \partial _{\nu }\overline{\omega }_{\mu }^{ac}\right) f^{abm}\left( D_{\nu }c\right) ^{b}\varphi _{\mu }^{mc}\right) \nonumber \\&+\int \mathrm {d}^d x\left( -\gamma ^{2}gf^{abc}A_{\mu }^{a}(\varphi _{\mu }^{bc}+\overline{\varphi }_{\mu }^{bc})\right. \nonumber \\&- \left. d\left( N^{2}-1\right) \gamma ^{4} \right) . \end{aligned}$$
(4)

The fields \(\varphi _{\mu }^{ac}\), \(\overline{\varphi }_{\mu }^{ac}\) are a set of bosonic fields each with \(d(N^2-1)^2\) components, the fields \(\omega _{\mu }^{ac}\), \(\overline{\omega }_{\mu }^{ac}\) are ghost fields. The Gribov mass parameter \(\gamma ^2\) is not free, but it is determined in a self-consistent way through the gap equation \(\left. \frac{\partial E}{\partial \gamma ^2}\right| _{\gamma ^2\ne 0}=0\), where \(E\) stands for the vacuum energy of the theory [57]. In more recent years, it was shown that nontrivial vacuum effects can further alter the action, leading to the so-called refined Gribov–Zwanziger (action) [811]. It is also worthwhile to mention that recently, progress has been made to generalize the Gribov–Zwanziger construction to the maximal Abelian gauge [1216]. Also studies of Gribov copies in curved space times were considered [17, 18], next to theories including Higgs scalar fields [19, 20] or supersymmetric gauge theories [21, 22].

An alternative way to cope with the Gribov ambiguity, based on sampling over all copies, has been put forward in [23]. The net result was the massive Curci–Ferrari action [24, 25], where the presence/value of the mass is, to our knowledge, a somewhat delicate issue related to taking a set of limits; see also [26] for similar questions. Nonetheless, a decent effective description of two- and higher point functions was feasible [2729].

In the Gribov–Zwanziger approach, we are implicitly assigning an equal weight to all further copies inside \(\Omega \) and sample over those. In the literature, this is sometimes referred to as the minimal Landau gauge. Such sampling might be ill-posed, if the number of copies inside \(\Omega \), \(\mathcal {N}[A]\), depends on the gauge field, but it is the best we can do analytically. Almost nothing is known about \(\mathcal {N}[A]\). Let us refer to [30] for a more elaborate discussion on the potential pitfalls of extending the Landau gauge beyond its perturbative definition.

It is perhaps also worthwhile to recall that \(\Omega \) itself is not free from residual gauge copies. This was discussed in mathematical detail in [3]; see also [31]. Ideally, one would have to resort to restricting the path integral to a region free of Gribov copies, defining an absolute Landau gauge. This so-called fundamental modular region (FMR) in the Landau gauge would correspond to the region \(\Lambda \subset \Omega \) built from the absolute (in contrast to the local) minima of the functional \(\mathcal {R}[A]\) (2). Unfortunately, there are no simple criteria which a generic \(A_\mu ^a\) should fulfill to belong to the FMR \(\Lambda \). Not much is known about \(\Lambda \), apart from the fact that it is also convex [3] and evidently, it is also bounded in every direction. Even numerically, it is extremely hard to locate (all) absolute minima of a functional in a complicated space. This is immediately the main reason why the Gribov–Zwanziger restriction to \(\Omega \) is mostly used in practical applications.

Even \(\Lambda \) itself is not free of copies, its boundary \(\partial \Lambda \) can display degenerate absolute minima as shown in [32]. It is also known that \(\Omega \) and \(\Lambda \) share part of their boundary [32]. In e.g. [33], it was conjectured and analytically motivated that \(\langle \ldots \rangle _\Omega =\langle \ldots \rangle _\Lambda \). Until now, this remains, however, a conjecture, and not all lattice studies comply with it [34]. In any case, interpreting lattice results with regard to the (continuum, thus infinite-volume) conjecture should be done with care, since the infinite volume is never reached.

The presence of copies inside the Gribov region \(\Omega \) opens the door for alternative completions of the Landau gauge; see also [35] for a discussion. For example, one could pick, per gauge configuration, as the relevant contribution to the path integral a copy strictly inside \(\Omega \) and not the one on the boundary \(\partial \Omega \) [36]. Gauge dependent quantities will in general be dependent on the specific copy selection.Footnote 1 Though, in principle, when evaluating gauge invariant quantities, it should be of no consequence whatever copy is used. Though, such selection procedure is at the practical level only feasible in the lattice context, and probably even only in at relatively small volumes, for which there is relative control on the set of Gribov copies.

In functional formalisms, with typical proponent the Dyson–Schwinger equations [3740], it was discussed in [41, 42] that the form/solutions of the equations should not change by cutting off the region of integration at the boundary \(\partial \Omega \).

In the following section, an analytical argument will be provided that—at least in the employed continuum setting—the path integral over the Gribov region \(\Omega \) is still dominated by the integration over its boundary \(\partial \Omega \) when making use of the Cucchieri–Mendes bound on the ghost propagator [43, 44]. In the third section, we focus on answering a question, originally posed by Maas [45, 46] and improved upon in [30], concerning the implementation at the functional level of a boundary condition on the (inverse) ghost propagator. Besides that, also some quantum aspects of the proposed action, including its renormalizability, are handled in the fourth section. We end with a more concise discussion of the obtained results.

2 The Cucchieri–Mendes bound on the ghost propagator and its implications

In [43], it was derived how the horizon function, defined via

$$\begin{aligned} h(A)&= \frac{1}{Vd(N^2-1)}g^2 \int \mathrm {d}^dx \mathrm {d}^dy f^{abc}A_\mu (x) \left[ \mathcal {M}^{-1}\right] ^{ad}\nonumber \\&\times f^{dec}A_\mu (y),\quad [V=\text {space-time volume}], \end{aligned}$$
(5)

obeys the following inequality:

$$\begin{aligned} h(A)\le \rho (A), \end{aligned}$$
(6)

where \(1-\rho (A)\) is the normalized distance of the gauge configuration \(A\) to \(\partial \Omega \), or, equivalently, \(\rho (A)\) is the normalized distance to \(0\in \Omega \). This distance can be defined on basis of the convexity of \(\Omega \) [1], in the sense that for each \(A\in \Omega \), we can consider in \(\Omega \) the ray connecting \(A_\mu \) with \(0\), calling \(B_\mu \) the intersection of that ray with \(\partial \Omega \). Doing so, \(\rho (A)\) corresponds to the scaling factor in the relation \(A_\mu = \rho (A) B_\mu \). In more mathematical terms [43], introducing the \(\mathcal {L}_2\) norm for vector valued SU(\(N\)) matrices, \(||A||^2=\frac{1}{dV(N^2-1)}\int \mathrm {d}^dx \text {Tr}\left[ A_\mu A_\mu ^\dagger \right] \), then \(\rho (A)= \frac{||A ||}{||B||}\). Notice that by construction \(\rho (A)\le 1\). Figure 1 displays the situation in a pictorial manner.

Fig. 1
figure 1

Pictorial representation of the Gribov region \(\Omega \) and the distance \(\rho (A)\)

In [47, 48], it was shown, order by order in a perturbative expansion, that the horizon function \(h(A)\) corresponds to the Gribov ghost form factor \(\sigma (k,A)\) at zero momentum. Imposing \(\sigma (0,A)\le 1\) [while keeping in mind that for the a posteriori constructed integration measure, \(\langle \sigma (k,A) \rangle \) is a decreasing function in \(k\)] is thus equivalent to imposing \(h(A)\le 1\). In the thermodynamic limit [5, 6], the inequality becomes an equality, see also [4], and using this we arrive at the horizon condition

$$\begin{aligned} \langle h(A) \rangle =1 \Leftrightarrow \langle \sigma (0,A) \rangle ^{1\mathrm{PI}}=1, \end{aligned}$$
(7)

where we added the superscript \(1\mathrm{PI}\) to make it clear that reference is made to the 1-particle-irreducible part of the ghost self-energy [48]. Evidently, the measure underpinning \(\langle \ldots \rangle \) in (7) is not the Faddeev–Popov one, since \(\sigma (0)\equiv \langle \sigma (0,A) \rangle ^{1\mathrm{PI}}\ne 1\) using standard perturbation theory. Let us also recall here that \(\sigma (0)=1\) leads to a more than quadratically IR divergent ghost propagator, a property first pointed out by Gribov [5], later on proven to all orders using a Ward identity argument by Zwanziger [7]. Naively, the Gribov–Zwanziger restriction boils down to taking into account only gauge configurations living on the Gribov horizon \(\partial \Omega \), where the ghost propagator diverges more strongly due to the presence of the zero modes.

Though, the observation of an IR enhanced ghost propagator is in disagreement with contemporary lattice simulationsFootnote 2 [4959], a fact later recognized in the Gribov–Zwanziger theory to be caused by a vacuum instability [811] that further changes the effective action in a dynamical fashion. Similar results were found using other approaches [27, 6063]. For a nonextensive list of recent additions to these research programs, let us refer to [6476].

Given the a priori sharper Cucchieri–Mendes bound (6), the possibility seems to be lurking around the corner for a geometrical reason behind a non-enhanced ghost in the infrared, given that one is actually integrating over a smaller subregion of the Gribov region \(\Omega \), since for many gauge configurations it will hold that \(\rho (A)<1\). Yet, depending on the number of configurations that do have \(\rho (A)=1\), ghost enhancement would still be allowed on mere classical/geometrical grounds but not realized explicitly because of another reason; in particular the dynamical refinement as discussed in [811], which would be a reason on quantum grounds. We refer to [43] for many more details, including lattice experiments and related discrete versions of the bound.

In the remainder of this section, we come to our first result: directly working in the continuum, we shall show that the path integral over the Gribov region \(\Omega \) is still dominated by the integration over the horizon \(\partial \Omega \). Indeed, consider the following path integral:

$$\begin{aligned}&\int [\mathcal {D} A] \theta [\rho (A)-h(A)] \theta [1-\rho (A)]\mathrm{e}^{-S_{\mathrm{FP}}[A]}\nonumber \\&\quad =\int _\Omega [\mathcal {D}A] \theta [\rho (A)-h(A)] \mathrm{e}^{-S_{\mathrm{FP}}[A]}=(*) \end{aligned}$$
(8)

which encompasses the Cucchieri–Mendes bound (6) via the first \(\theta \)-function, next to the restriction to the Gribov region via the second \(\theta \)-function. Now given that the region \(\Omega \) is for sure star-convex, it is natural to use “polar coordinates” to rewrite the integral as an integral over the boundary \(\Omega \) times a “normalized radial” integral;Footnote 3

$$\begin{aligned} (*)= \int _{\partial \Omega } [\mathcal {D}B] \int _0^{1} \mathrm {d}\rho \rho ^N \theta [\rho -h(\rho B)] \mathrm{e}^{-S_{\mathrm{FP}}[\rho B]}, \end{aligned}$$
(9)

where \(N=\dim (\Omega )-1=\infty \). We first rewrite the \(B\)-integral via the identification

$$\begin{aligned} \int _{\partial \Omega } [\mathcal {D}B](\ldots )=\int [\mathcal {D}B] \delta (h(B)-1)(\ldots ). \end{aligned}$$
(10)

The latter integral can be interpreted as defining a microcanonical average. Making use of the thermodynamic limit [6, 77], the former averages equal those generated by a Boltzmann distribution (see also later for some more details), defined by

$$\begin{aligned} \int _{\partial \Omega } [\mathcal {D}B](\ldots )&= \int [\mathcal {D}B] \delta (h(B)-1)(\ldots )\nonumber \\&\rightarrow \int [\mathcal {D} B] \mathrm{e}^{\gamma ^4\int \mathrm {d}^4x h(B)}(\ldots ), \end{aligned}$$
(11)

with \(\gamma ^4\) the critical (Gribov) parameter determined by the self-consistency equation

$$\begin{aligned} \langle h \rangle =1, \end{aligned}$$
(12)

where the measure adopted is the one corresponding to the r.h.s. of (11). We thus get for the partition function

$$\begin{aligned} (*)\!=\! \int [\mathcal {D} B] \int _0^{1} \mathrm {d}\rho \rho ^N \theta [\rho -h(\rho B)] \mathrm{e}^{-S_{\mathrm{FP}}[\rho B]\!+\!\gamma ^4\int \mathrm {d}^4x h(B)}.\nonumber \\ \end{aligned}$$
(13)

Then we use the Laplace approximation rule for an integral of the form

$$\begin{aligned} \lim _{N\rightarrow \infty } \int _0^1 \mathrm {d}xx^N f(x)&= \lim _{N\rightarrow \infty } \int _0^1 \mathrm {d}x \mathrm{e}^{N\ln x} f(x)\nonumber \\&\propto \lim _{N\rightarrow \infty } \mathrm{e}^{N\ln x_0} f(x_0), \end{aligned}$$
(14)

where \(x_0=1\) is the maximum of \(\ln x\) over the interval \([0,1]\). Thus, it is clear the \(\rho \)-integral in (13) will be reduced to its value at \(\rho =1\) (i.e. on the horizon). We find for our partition function

$$\begin{aligned} (*)= \int [\mathcal {D} B] \theta [1-h(B)] \mathrm{e}^{-S_{\mathrm{FP}}[B]+\gamma ^4\int \mathrm {d}^4x h(B)}. \end{aligned}$$
(15)

In the latter expression, we recognize almost the standard Gribov–Zwanziger action. The difference lies in the presence of the \(\theta [1-h(B)]\), but this is actually an obsolete factor. Indeed, thanks to the thermodynamic limit/infinite-dimensionality of the configuration space [4] the \(\theta \)-function collapses into a \(\delta \)-function,Footnote 4 and we can repeat the microcanonical \(\rightarrow \) Boltzmann averaging. Though, since the measure \([\mathrm {d}B] \mathrm{e}^{-S_{\mathrm{FP}}[B]+\gamma ^4\int d^4x h(B)}\) already ensures \(\langle h \rangle =1\), there is no need for an additional critical parameter \(\hat{\gamma }^4\) (more precisely, \(\hat{\gamma }^4=0\)), and we eventually arrive at

$$\begin{aligned} (*)= \int [\mathcal {D}B] \mathrm{e}^{-S_{\mathrm{FP}}[B]+\gamma ^4\int \mathrm {d}^4x h(B)} \end{aligned}$$
(16)

which is nothing else than the Gribov–Zwanziger action in its nonlocal formulation.

It should not have come as a surprise that imposing the somewhat more stringent Cucchieri–Mendes bound leads to exactly the Gribov–Zwanziger action after all. In the original derivation, the argument is that the dominant integration region is the boundary of \(\Omega \). Naively speaking, imposing the Cucchieri–Mendes bound via (8), amounts to integrating over an outer shell of the ellipsoid \(\Omega \). The dominant region will still be the outer boundary of this infinite-dimensional shell, being \(\partial \Omega \).

Let us here also point to [7880] for a numerical verification that the relevant configurations are located close to \(\partial \Omega \), as it was found that the average value of the smallest nontrivial eigenvalue of the Faddeev–Popov operator approaches zero for growing lattice volume.

3 Using the inverse ghost form factor as a “nonperturbative gauge parameter”: construction of the corresponding action

In this section, we will take a closer look at the so-called Landau \(B\)-gauges as first introduced in [45, 46] and studied afterwards in [59], albeit mostly in a numerical lattice context. Obviously, the Landau gauge is assumed to begin with. The motivation lies in the Dyson–Schwinger quantum equations of motion (DSEs). These functional equations appear to have multiple solutions, and finding a unique solution depends on a choice of boundary conditions; see e.g. [81, 82]. In the context of the DSEs, the (not always explicitly mentioned) boundary condition that is imposed is making a choice for the \(1\mathrm{PI}\) ghost self-energy at zero momentum; see e.g. [83, Sect V.C]. This observation was explicitly made and numerically tested in [84] for the ghost propagator in the Coulomb gauge.

More precisely, in the Landau gauge, using global color invariance, one may write

$$\begin{aligned} \mathcal {G}^{ab}(p^2)=\frac{g(p^2)}{p^2}\delta ^{ab}. \end{aligned}$$
(17)

One then imposes the requirement that the ghost form factor \(g(p^2)\) obey

$$\begin{aligned} g^{-1}(0)=B. \end{aligned}$$
(18)

A special case is provided by \(B=0\), which corresponds to a ghost propagator that is more singular than \(1/p^2\) in the deep infrared. This is commonly known as the “scaling” scenario, while \(B>0\) corresponds to the “decoupling/massive” scenario, the latter being realized in current lattice simulations.Footnote 5 This terminology was coined in [62]. Although the question whether both scenarios could be realized on the lattice is to our knowledge not yet answered with strict rigor, various reasons why finding a scaling solution on the lattice would be rather surprising can be found in the discussion of [43]; see also [85] for an argumentation in the continuum. Interestingly, the latter paper also discusses the effective action related to both classes of solutions, giving evidence that the massive/decoupling solution corresponds to a more stable situation. The fact that multiple distinct solutions do exist, but that the theory itself singles one out on energetic grounds (smaller vacuum energy) sounds, evidently, physically appealing and logical. Arguments suggesting the possible existence of the scaling solution were provided in [86].

Our goal now is to show that it is possible to enforce the specific boundary condition \(g^{-1}(0)=B\) directly at the level of the action, without making compromises w.r.t. the renormalizability of the construction. This gives an explicit answer to the speculations of [45, 46] if it is possible to implement at the functional level the boundary condition. Notice that being able to implement all constraints at the level of the Lagrangian/action allows for a clean discussion of the properties of the eventual solution. We will therefore combine the tools of [48, 77]. Let us consider the Faddeev–Popov operator, \(\mathcal {M}^{ab}\) and for the moment we treat \(A_\mu \) as a classical external field. The ghost propagator in Fourier space reads

$$\begin{aligned} \langle \bar{c}^a(p) c^b(-q) \rangle&= \int \mathrm {d}^d x \; \int \mathrm {d}^d y \langle \bar{c}^a(x) c^b(y) \rangle \mathrm{e}^{-ipx}\mathrm{e}^{iqy},\nonumber \\ \end{aligned}$$
(19)

with

$$\begin{aligned} \langle \bar{c}^a(x) c^b(y) \rangle&= \frac{\int [\mathcal{D}c][\mathcal{D}\bar{c}] \;\bar{c}^a(x) c^b(y) \mathrm{e}^{\int \bar{c}^c \mathcal{M}^{cd} c^d }}{\int [\mathcal{D}c][ \mathcal{D} \bar{c}]\;\mathrm{e}^{\int \bar{c}^c \mathcal{M}^{cd} c^d}}. \end{aligned}$$
(20)

For the ghost form factor, a compact expression can be derived if we look at the color trace of (20), that is,

$$\begin{aligned} \mathcal {G} (p,A)&= \frac{1}{V(N^2-1)} \langle \bar{c}^a(p) c^a(-q) \rangle \delta (p+q)\nonumber \\&= \frac{1}{p^2} (1 + \sigma (p,A)), \end{aligned}$$
(21)

with \(V\equiv \delta (0)\) the (infinite) 4d spacetime volume. Using Wick’s theorem, the transversality of the Landau gauge, and some algebra, it was shown in [48] that \(\sigma (p=0,A)\) can be compactly rewritten as

$$\begin{aligned} \sigma (0, A)&= -\frac{g^2}{VD(N^2-1)} \int \frac{\mathrm {d}^D p}{(2\pi )^D} \nonumber \\&\times \int \frac{\mathrm {d}^D q}{(2\pi )^D} A^{ab}_{\mu }(-p) \left( \mathcal{M}^{-1}\right) ^{bc}(p,q)A^{ca}_{\mu }(q)\nonumber \\&\equiv \frac{H(A)}{DV(N^2-1)}, \end{aligned}$$
(22)

where \(A_\mu ^{ab}=-f^{abc}A_\mu ^c\) and \(\mathcal {M}^{ab}(p,q)= q^2\delta ^{ab}\delta (p-q) - gA^{ab}_{\mu }(p-q)iq_{\mu }\).

Now, we can quantize the gauge field as well, and (22) can be translated into the following equality:

$$\begin{aligned} \mathcal{G} (k)&= \langle \mathcal{G} (k,A)\rangle ^{\mathrm{conn}} = \frac{1}{k^2} (1 + \langle \sigma (k,A)\rangle ^{\mathrm{conn}}) \nonumber \\&= \frac{1}{k^2} \frac{1}{(1 - \langle \sigma (k,A)\rangle ^{1\mathrm{PI}})}, \end{aligned}$$
(23)

where “conn” stands for the connected set of diagrams and \(1\mathrm{PI}\) for the 1-particle-irreducible ones. The boundary condition (18) is then equivalent to requiring

$$\begin{aligned} \langle \sigma (0,A) \rangle ^{1\mathrm{PI}}=1-B\equiv \hat{B}. \end{aligned}$$
(24)

How to impose this at the level of the partition function? This can be accommodated for by a replacement of the Landau gauge Faddeev–Popov measure as follows:

$$\begin{aligned} \int [\mathcal {D}\mu _{\mathrm{FP}}]&\equiv \int [\mathcal {D}A] \det \mathcal {M}\delta (\partial A)\mathrm{e}^{-S_{\mathrm{YM}}}\rightarrow \int [\mathcal {D}\mu ']\nonumber \\&\equiv \int [\mathcal {D}\mu _{\mathrm{FP}}]\delta (\sigma (0)-\hat{B}) \nonumber \\&= \int [\mathcal {D}A] \det \mathcal {M}\delta (\partial A)\delta (H(A)\nonumber \\&-\hat{B} dV(N^2-1))\mathrm{e}^{-S_{\mathrm{YM}}}. \end{aligned}$$
(25)

In statistical mechanics terms, the latter measure corresponds to taking a microcanonical ensemble average. In the thermodynamic limit, \(\sharp \) d.o.f.\(~\rightarrow \infty \), \(V\rightarrow \infty \); \(\sharp \) d.o.f./\(V\) = constant; the measure can be reformulated with a steepest descent evaluation (a.k.a. saddle point approximation) to arrive at a Boltzmann ensemble.Footnote 6 More precisely [5, 77],

$$\begin{aligned}&\int [\mathcal {D}\mu _{\mathrm{FP}}]\delta (H(A)-\hat{B} dV(N^2-1))\nonumber \\&\quad =\int [\mathcal {D}\mu _{\mathrm{FP}}] \int _{-i\infty +\epsilon }^{i\infty +\epsilon }\frac{\mathrm {d}\beta }{2\pi i}\mathrm{e}^{-\beta (H(A)-\hat{B} dV(N^2-1))}\nonumber \\&\quad =\int _{-i\infty +\epsilon }^{i\infty +\epsilon }\frac{\mathrm {d}\beta }{2\pi i} \mathrm{e}^{-E(\beta )}, \end{aligned}$$
(26)

with

$$\begin{aligned} E(\beta )\!=\!-\ln f(\beta )\!=\!-\ln \int [\mathcal {D}\mu _{\mathrm{FP}}]\mathrm{e}^{-\beta (H(A)-\hat{B} dV(N^2-1))}.\nonumber \\ \end{aligned}$$
(27)

Power counting shows that \(\beta \) has mass dimension \(4\). To evaluate the \(\beta \)-integral of (26), we need to solve

$$\begin{aligned}&\left. \frac{\partial E(\beta )}{\partial \beta }\right| _{\beta =\beta ^*}= \frac{f'(\beta ^*)}{f(\beta ^*)}=0 \end{aligned}$$
(28)

to locate the critical point \(\beta ^*\) and make sure we can deform the contour to pass through \(\beta ^*\). This is allowed if the integration along the imaginary axis can be extended to a contour integration. From the representation (26) it is actually clear that \(E(\beta )\) corresponds to the vacuum energy up to the (infinite) volume factor \(V\). On dimensional grounds, in \(d\) dimensions, we will have \(E(\beta )\propto \beta \), so if we close the contour on the right (\(Re(\beta )>0\)), the \(\int \mathrm{e}^{-E(\beta )}\) will vanish for \(\beta \rightarrow \infty \) and we can safely deform the contour to pass through \(\beta ^*\), as long as \(\beta ^*\ge 0\) and if there are no poles or cuts located in the complex \(\beta \)-plane with \(Re(\beta )>0\). We will come back to the latter issue further on. A little algebra results in \(\frac{\partial ^2 E}{\partial \beta ^2}=\langle (H-\langle H \rangle )^2 \rangle >0\), so we expect that \(\frac{\partial ^2 E}{\partial \beta ^2}>0\) at the critical point. This will correspond to a global minimum of the vacuum energy,Footnote 7 and the steepest descent evaluation should thus be allowed.

Returning to the saddle point equation (28), we notice that it can be rephrased as

$$\begin{aligned} \langle H(A) \rangle =\hat{B} dV(N^2-1), \end{aligned}$$
(29)

where the expectation value is taken w.r.t. the new measure

$$\begin{aligned} \int [\mathcal {D}A] \det \mathcal {M}\delta (\partial A)\mathrm{e}^{-S_{\mathrm{YM}}}\mathrm{e}^{-\beta ^*(H(A)-\hat{B} dV(N^2-1))}. \end{aligned}$$
(30)

Now, since \(E\) is composed of vacuum bubble contributions, any potential one-particle reducible diagram contributing to \(E\) must be proportional to the square of the zero-momentum expectation value of the field propagating in the reducible line. This refers to nothing else than the vacuum expectation value of the fundamental field. Due to Lorentz and/or global color invariance, such condensates are, however, zero. Consequently, the relation (29) actually reads

$$\begin{aligned} \langle H(A) \rangle ^{1\mathrm{PI}}=\hat{B} dV(N^2-1), \end{aligned}$$
(31)

which, upon using (22) again, leads to

$$\begin{aligned} \langle \sigma (0,A) \rangle ^{1\mathrm{PI}}=\hat{B}, \end{aligned}$$
(32)

that is, we indeed have constructed a novel partition function that precisely gives the desired result (18), \(g^{-1}(0)=B\).

A word of caution is in order here. Writing down (25), we made the a priori assumption that for each gauge configuration, there is at least one gauge equivalent configuration leading to the imposed value of \(B\) while satisfying the Landau gauge condition.Footnote 8 Unfortunately, there is no formal proof of this. To our knowledge, at most it is rigorously known that each gauge orbit has at least one representant inside the Gribov region \(\Omega \), including its boundary [2]. We will, however, stick to the assumption that every value of \(B\ge 0\) is allowed. This is in part motivated by the (less strong) assumption of the Dyson–Schwinger paper [62] that, after taking expectation values, the ghost propagator’s form factor can attain any \(B>0\) as its zero momentum value. Further numerical studies concerning the range of allowed \(B\)’s and related discussion can be found in [30, 35, 43, 59, 87].

To bring the partition function in a more usual, i.e. exponential, form, we can, completely analogously to the original Zwanziger derivation, introduce a set of auxiliary bosonic \(\varphi _\mu ^{ab}\), \(\overline{\varphi }_\mu ^{ab}\), and fermionic fields, \(\omega _\mu ^{ab}\), \(\overline{\omega }_\mu ^{ab}\), to localize the nonlocal operator \(H(A)\). We eventually find that the partition function is completely determined in terms of the action

$$\begin{aligned} S&= \frac{1}{4}\int \mathrm {d}^d x\; F^a_{\mu \nu } F^a_{\mu \nu }+\int \mathrm {d}^d x\,\left( b^a \partial _\mu A_\mu ^a +\overline{c}^a \partial _\mu D_\mu ^{ab} c^b \right) \nonumber \\&+\int \mathrm {d}^d x \left( \overline{\varphi }_{\mu }^{ac} \partial _{\nu } D_\nu ^{am}\varphi _{\mu }^{mc} -\overline{\omega }_{\mu }^{ac} \partial _{\nu } D_\nu ^{am} \omega _{\mu }^{mc} \right. \nonumber \\&-\left. g\left( \partial _{\nu }\overline{\omega }_{\mu }^{ac}\right) f^{abm}\left( D_{\nu }c\right) ^{b}\varphi _{\mu }^{mc}\right) \nonumber \\&+\int \mathrm {d}^d x\left( -\sqrt{\beta ^*}gf^{abc}A_{\mu }^{a}\left( \varphi _{\mu }^{bc}+\overline{\varphi }_{\mu }^{bc}\right) \right. \nonumber \\&\left. - d\hat{B}\left( N^{2}-1\right) \beta ^* \right) . \end{aligned}$$
(33)

Indeed, being Gaussian, both the \(\varphi \)- and the \(\omega \)-integrals can be executed exactly, their respective determinants canceling due to their opposite Grassmann parity.

The mass scale \(\beta ^*\) is still to be fixed from the gap equation (28), which can, in the local formulation, also be expressed as

$$\begin{aligned} \left\langle gf^{abc}A_{\mu }^{a}\left( \varphi _{\mu }^{bc}+\overline{\varphi }_{\mu }^{bc}\right) \right\rangle =2d\hat{B}(N^2-1)\sqrt{\beta ^*}. \end{aligned}$$
(34)

Obviously, the action we just constructed is a generalization of the Gribov–Zwanziger one, which is recovered in the limit \(\hat{B}=1\), in which case the original no-pole condition of Gribov is implemented.

From the local formulation, to which we can apply the usual effective action formalism, it becomes evident that the vacuum energy, upon including the quantum corrections, will be a function of the type \(E\left( \beta ^*, g^2,\ln \frac{\sqrt{\beta ^*}}{\mu ^2}\right) \). Indeed, apart from the pure vacuum term, \(\sqrt{\beta ^*}\) only appears in the quadratic term that mixes the \(A_\mu ^a\)-field with the \(\varphi _\mu ^{ac}\), \(\overline{\varphi }_\mu ^{ac}\) fields, thus it will only enter the (mixed) \((b^a,A_\mu ^a, \varphi _\mu ^{ac}, \overline{\varphi }_\mu ^{ac})\) propagators. In the absence of any other scale, up to the global prefactor, the only way it can enter the vacuum energy is logarithmically. Taking as before \(Re(\beta )>0\), we avoid the cut of the logarithm, and we can safely deform the contour for the saddle point evaluation.

By a simple rescaling \(\beta ^*\rightarrow \frac{\beta ^*}{\hat{B}}\), the “gauge parameter” \(\hat{B}\) will only enter via the \((b^a, A_\mu ^a, \varphi _\mu ^{ac}, \overline{\varphi }_\mu ^{ac})\) propagators, a situation rather akin to one with a genuine gauge parameter. Though the situation is very different, in particular physical quantities, as \(E\), may explicitly depend on \(\hat{B}\); see later in this paper.

4 The inverse ghost form factor as a “nonperturbative gauge parameter”: some quantum aspects

4.1 Renormalizability

In order to discuss the renormalizability, we can follow the same argumentation as for the Gribov–Zwanziger action, as the only difference is the tree level pure vacuum term, which leads to a different gap equation. We use the setup and notations ofFootnote 9 [88] which fits within the general algebraic formalism [91]. We set \(\sqrt{\beta }=\gamma ^2\) and first introduce a multi-index \(i\equiv (\mu ,a)\) which allows one to rewrite the action (33) as

$$\begin{aligned} S&= S_{0} + S_{\gamma } , \end{aligned}$$
(35)

with

$$\begin{aligned} S_{0}&= S_\mathrm {YM}+ S_\mathrm {gf}+ \int \mathrm {d}^d x \left( \overline{\varphi }_i^a \partial _\mu \left( D_\mu ^{ab} \varphi ^b_i \right) \right. \nonumber \\&-\left. \overline{\omega }_i^a \partial _\mu \left( D_\mu ^{ab} \omega _i^b \right) - g f^{abc} \partial _\mu \overline{\omega }_i^a D_\mu ^{bd} c^d \varphi _i^c \right) \nonumber , \\ S_{\gamma }&= -\gamma ^{2}g\int \mathrm {d}^{d}x\left( f^{abc}A_{\mu }^{a}\varphi _{\mu }^{bc} +f^{abc}A_{\mu }^{a}\overline{\varphi }_{\mu }^{bc} \right. \nonumber \\&+ \left. \hat{B}\frac{d}{g}\left( N^{2}-1\right) \gamma ^{2} \right) . \end{aligned}$$
(36)

In order to find a sufficiently large set of Ward identities, it is useful to embed the previous action into another, more general one,

$$\begin{aligned} \tilde{\Sigma }&= S_0 + S_\sigma , \end{aligned}$$
(37)

whereby

$$\begin{aligned} S_\sigma&= s\int \mathrm {d}^d x \left( -U_\mu ^{ai} D_\mu ^{ab} \varphi _i^b - V_\mu ^{ai} D_{\mu }^{ab} \overline{\omega }_i^{b} - U_\mu ^{ai} V_\mu ^{ai} \right. \nonumber \\&+\left. T_\mu ^{a i} g f_{abc} D^{bd}_\mu c^d \overline{\omega }^c_i \right) \nonumber \\&= \int \mathrm {d}^d x \left( -M_\mu ^{ai} D_\mu ^{ab} \varphi _i^b - gf^{abc} U_\mu ^{ai} D^{bd}_\mu c^d \varphi _i^c \right. \nonumber \\&+\left. U_\mu ^{ai} D_\mu ^{ab} \omega _i^b - N_\mu ^{ai} D_\mu ^{ab} \overline{\omega }_i^b - V_\mu ^{ai} D_\mu ^{ab} \overline{\varphi }_i^b \right. \nonumber \\&+\left. gf^{abc} V_\mu ^{ai} D_\mu ^{bd} c^d \overline{\omega }_i^c -\hat{B} M_\mu ^{ai} V_\mu ^{ai}+\hat{B} U_\mu ^{ai} N_\mu ^{ai} \right. \nonumber \\&+\left. R_\mu ^{ai} g f^{abc} D_\mu ^{bd} c^d \overline{\omega }^c_i + T_\mu ^{ai} g f_{abc} D^{bd}_\mu c^d \overline{\varphi }^c_i\right) . \end{aligned}$$
(38)

The Yang–Mills fields are subject to their standard BRST variation,

$$\begin{aligned} s A_\mu ^a&= -D_{\mu }^{ab}c^b, \quad s c^a = \frac{g}{2}f^{abc}c^b c^c, \nonumber \\ s \overline{c}^a&= b^a, \quad s b^a =0, \end{aligned}$$
(39)

while the auxiliary fields are BRST doublets,

$$\begin{aligned} s \varphi _i^a&= \omega _i^a, \quad s \omega _i^a = 0, \nonumber \\ s \overline{\omega }_i^a&= \overline{\varphi }_i^a, \quad s \overline{\varphi }_i^a =0. \end{aligned}$$
(40)

We have also introduced three new BRST doublets of sources (\(U_\mu ^{ai}\), \(M_\mu ^{ai}\)), (\(V_\mu ^{ai}\), \(N_\mu ^{ai}\)), and (\(T_\mu ^{ai}\), \(R_\mu ^{ai}\)) transforming according to

$$\begin{aligned} sU_{\mu }^{ai}&= M_{\mu }^{ai}, \quad sM_{\mu }^{ai}=0, \nonumber \\ sV_{\mu }^{ai}&= N_{\mu }^{ai}, \quad sN_{\mu }^{ai}=0,\nonumber \\ sT_{\mu }^{ai}&= R_{\mu }^{ai}, \quad sR_{\mu }^{ai}=0. \end{aligned}$$
(41)

For the special set of source values

$$\begin{aligned}&U_\mu ^{ai} = N_\mu ^{ai} = T_\mu ^{ai} = 0 , \nonumber \\&M_{\mu \nu }^{ab}= V_{\mu \nu }^{ab}= R_{\mu \nu }^{ab} = \gamma ^2 \delta ^{ab}\delta _{\mu \nu } , \end{aligned}$$
(42)

we recover the action of interest, \(S\) (35). To consistently define the nonlinear BRST variations of \(A_\mu ^a\) and \(c^a\), we also temporarily need to address

$$\begin{aligned} \Sigma&= S_\sigma + \int \mathrm{d}^4x \left( -K_\mu ^a D_\mu ^{ab}c^b+ \frac{g}{2}L^a f^{abc}c^b c^c\right) , \end{aligned}$$
(43)

with

$$\begin{aligned} sK_{\mu }^{a}&= 0, \quad sL^a=0. \end{aligned}$$
(44)

The “full” action \(\Sigma \) obeys several powerful Ward identities. We refrain from writing them all down there, they can be found in the Appendices of [88]; the only difference with the action discussed there is the presence of the additional parameter \(\hat{B}\), which only affects the pure source (vacuum) term \(\propto ~MV-UN\). Clearly, a change of the latter will never have any consequence whatsoever on the UV structure of the theory, so the theory for any value of \(\hat{B}\) will be renormalizable, given that the renormalizability has already been established in e.g. [7, 88, 92] for \(\hat{B}=1\). In fact, there is an additional Ward identity directly related to \(\hat{B}\), namely

$$\begin{aligned} \frac{\partial \Sigma }{\partial \hat{B}}&= \int \mathrm {d}^4x\left( -M_\mu ^{ai}V_\mu ^{ai} +U_\mu ^{ai}N_\mu ^{ai}\right) . \end{aligned}$$
(45)

The latter is a pure source term, as such the identity remains valid at the quantum level for the \(1\mathrm{PI}\) effective action \(\Gamma \). Consequently, as \(\hat{B}\) does appear in the bare action but not in the counterterm, the renormalization of this pure source term, \(\hat{B}(-MV + UN)\), will happen as

$$\begin{aligned}&Z_{\hat{B}}\left( Z_{M}Z_{V}\hat{B}MV - Z_{U}Z_{N}\hat{B}UN\right) \nonumber \\&\quad = \left( 1-h \frac{a_{1}}{\hat{B}}\right) \hat{B}(-MV + UN), \end{aligned}$$
(46)

where \(h\) is the infinitesimal parameter associated to the quantum perturbation of the classical action [91]. As can be inferred from [88], \(Z_{M}Z_{V}=Z_{U}Z_{N}=(1- h a_{1})=Z^{2}_{\gamma ^{2}}\) and so we have

$$\begin{aligned} Z_{\hat{B}} = \left[ 1 - h a_{1} \left( \frac{1}{\hat{B}}-1\right) \right] , \end{aligned}$$
(47)

which finally gives

$$\begin{aligned} Z_{\hat{B}} = \left( 1 - h a_{1} \right) ^{\left( \frac{1}{\hat{B}}-1\right) } = Z^{2\left( \frac{1}{\hat{B}}-1\right) }_{\gamma ^{2}} . \end{aligned}$$
(48)

In the above, \(a_{1}\) is an arbitrary parameter that comes from the counterterm and \(Z_{\gamma ^{2}}\), which is the multiplicative renormalization factor for \(\gamma ^2\), is given by \(Z^{2}_{\gamma ^{2}} = \left( 1+\frac{3}{2}\frac{g^{2}N}{16\pi ^{2}}\frac{1}{\varepsilon }\right) \) as can be learned from [9] (notations and conventions are as follows [88]).

The only possible change in the other Ward identities of [88] will at most occur in the pure source term, that is, in the potential harmless linear breaking terms (even zeroth order terms) of the classical/quantum Ward identities. For further usage, we will need the following Ward identities:

$$\begin{aligned}&\frac{\delta \Sigma }{\delta b^a}= \partial _\mu A_\mu ^a ,\end{aligned}$$
(49)
$$\begin{aligned}&\frac{\delta \Sigma }{\delta \omega _i^a}+\partial _\mu \frac{\delta \Sigma }{\delta N_\mu ^{ai}}-gf^{abc}\overline{\omega }_i^b\frac{\delta \Sigma }{\delta b^c}= -\hat{B} \partial _\mu U_\mu ^{ai}\nonumber \\&\quad +\,D_\mu ^{ab}U_\mu ^{ib} \end{aligned}$$
(50)

andFootnote 10

$$\begin{aligned}&\frac{\delta \Sigma }{\delta \overline{\omega }_i^a}+\partial _\mu \frac{\delta \Sigma }{\delta U_\mu ^{ia}}+gf^{abc}\left( V_\mu ^{ci}+R_\mu ^{ci}\right) \frac{\delta \Sigma }{\delta K_\mu ^b}\nonumber \\&\quad = \hat{B} \partial _\mu N_\mu ^{ai}-D_\mu ^{ab}N_\mu ^{bi}, \end{aligned}$$
(51)

next to the global invariance (thanks to the multi-index)

$$\begin{aligned} U_{ij}\Sigma&= 0\nonumber \\ U_{ij}&= \int \mathrm {d}^dx\left( \varphi _{i}^{a}\frac{\delta }{\delta \varphi _{j}^{a}}-\overline{\varphi }_{j}^{a}\frac{\delta }{\delta \overline{\varphi }_{i}^{a}}+\omega _{i}^{a}\frac{\delta }{\delta \omega _{j}^{a}}-\overline{\omega }_{j}^{a}\frac{\delta }{\delta \overline{\omega }_{i}^{a}} \right. \nonumber \\&- M^{aj}_{\mu } \frac{\delta }{\delta M^{ai}_{\mu }} -U^{aj}_{\mu }\frac{\delta }{\delta U^{ai}_{\mu }} + N^{ai}_{\mu }\frac{\delta }{\delta N^{aj}_{\mu }}\nonumber \\&+\left. V^{ai}_{\mu }\frac{\delta }{\delta V^{aj}_{\mu }} + R^{aj}_{\mu }\frac{\delta }{\delta R^{ai}_{\mu }} + T^{aj}_{\mu }\frac{\delta }{\delta T^{ai}_{\mu }} \right) , \end{aligned}$$
(52)

by means of which we can introduce the “\(i\)-charge” \(Q=U_{ii}\). Last but not least, we will also need the Slavnov–Taylor identity:

$$\begin{aligned}&\int d^4x\left( \frac{\delta \Sigma }{\delta K_\mu ^a}\frac{\delta \Sigma }{\delta A_\mu ^a}+\frac{\delta \Sigma }{\delta L^a}\frac{\delta \Sigma }{\delta c^a}+b^a\frac{\delta \Sigma }{\delta \overline{c}^a}+\overline{\varphi }_i^a\frac{\delta \Sigma }{\delta \overline{\omega }_i^a}\right. \nonumber \\&\quad \quad +\left. \omega _i^a\frac{\delta \Sigma }{\delta \varphi _i^a}+M_\mu ^{ai}\frac{\delta \Sigma }{\delta U_\mu ^{ai}}+N_\mu ^{ai}\frac{\delta \Sigma }{\delta V_\mu ^{ai}}+R_\mu ^{ai}\frac{\delta \Sigma }{\delta T_\mu ^{ai}}\right) \nonumber \\&\quad = 0. \end{aligned}$$
(53)

4.2 The anomalous dimension \(\gamma _{\hat{B}}\)

It is interesting to take a closer look at the anomalous dimension of \(\hat{B}\), which we define via

$$\begin{aligned} \gamma _{\hat{B}}= - \frac{\partial \ln Z_{\hat{B}}}{\partial \ln \bar{\mu }} \end{aligned}$$
(54)

and from (48) we learn

$$\begin{aligned} \ln Z_{\hat{B}} = \left( \frac{1}{\hat{B}}-1\right) \frac{3}{2}\frac{\alpha }{\varepsilon }, \quad \text {with } \alpha = \frac{g^{2}N}{(4\pi )^{2}}. \end{aligned}$$
(55)

Therefore we can write

$$\begin{aligned} \gamma _{\hat{B}}= - \frac{\partial \ln Z_{\hat{B}}}{\partial \alpha } \frac{\partial \alpha }{\partial \ln \bar{\mu }}. \end{aligned}$$
(56)

Setting \(g_{0}^{2} = \bar{\mu }^{\varepsilon }Z^{2}_{g}g^{2}\) we have

$$\begin{aligned} \frac{\partial \alpha }{\partial \ln \bar{\mu }} = \left( -\varepsilon g^{2}+2\beta (g)\right) \frac{N}{(4\pi )^{2}}, \end{aligned}$$
(57)

so that

$$\begin{aligned} \gamma _{\hat{B}}= - \left( -\varepsilon g^{2}+2\beta (g)\right) \frac{N}{(4\pi )^{2}} \left( \frac{1}{\hat{B}}-1\right) \frac{3}{2}\frac{1}{\varepsilon }, \end{aligned}$$
(58)

leading to

$$\begin{aligned} \gamma _{\hat{B}}= \left( \frac{1}{\hat{B}}-1\right) \frac{3}{2}\frac{g^{2}N}{(4\pi )^{2}}. \end{aligned}$$
(59)

According to (59) and since \(\hat{B} \ge 1\) the anomalous dimension of \(\hat{B}\) is thus positive and it has a fixed point at \(\hat{B}=1\), reached in the UV regime. The fact that \(\hat{B}=1\) is a fixed point to all ordersFootnote 11 is a cross-check that the original Gribov–Zwanziger action is a stable special case of the general action that gives a prescribed zero-momentum ghost form factor. Since \(\hat{B}=1\) is reached in the ultraviolet regime, it would appear that, starting from a choice of \(\hat{B}\) that does not correspond to the original Gribov–Zwanziger case (scaling solution), one will likely not end up with \(\hat{B}=1\) at low momentum.

4.3 Explicit proof of the boundary condition being realized

In this subsection, we provide a proof that the ghost propagator indeed obeys the boundary condition (18) when quantum corrections are taken into account. The argumentation will solely rest on the underlying Ward identities defining the quantum effective action. In the case of \(\hat{B}=1\), a proof can already be found in [7], but the more general new proof here presented is shorter and somewhat less technical and applies to general \(\hat{B}\). The gain is coming from explicit use of the quantum Slavnov–Taylor identity.

We depart from the partition function associated with the action (33). By adding sources coupled to the Faddeev–Popov ghosts and Gribov–Zwanziger ghosts, integrating them out exactly as they appear at most quadratically, one finds the general identification (see e.g. Section III of [88], albeit the proportionality factor mentioned there is now corrected here).

$$\begin{aligned}&\langle f(A) (\mathcal {M}^{-1})^{ab}(x,y) \rangle = \langle f(A) c^a(x)\overline{c}^b(y) \rangle \nonumber \\&\quad =-\frac{1}{4(N^2-1)} \left\langle f(A) \omega _i^a(x) \overline{\omega }_i^b(y)\right\rangle , \end{aligned}$$
(60)

where \(f(A)\) is a generic functional of the gauge fields; the \(4(N^2-1)\) corresponds to \(\delta _{ii}\), where \(i\) runs over \(d(N^2-1)\) values, with \(d\) the spacetime dimensionality. Rather than considering the ghost form factor at zero momentum itself, we can thus equally well look at the equivalent expression

$$\begin{aligned} g(0)&= -\frac{1}{4(N^2-1)^2} \lim _{p^2\rightarrow 0} p^2 \int \frac{\mathrm {d}^4x}{(2\pi )^4} \mathrm{e}^{ip(x-y)} \nonumber \\&\times \left\langle \omega _i^a(x) \overline{\omega }_i^a(y)\right\rangle . \end{aligned}$$
(61)

For now, let us focus on the Green function \(\langle \omega _i^a(x) \overline{\omega }_i^a(y) \rangle \), or even more specifically, to its \(1\mathrm{PI}\) counterpart (inverse propagator). The latter can be inferred from

$$\begin{aligned} \left. \frac{\delta ^2 \Gamma }{\delta \omega _i^a(x) \delta \overline{\omega }_j^c(y)}\right| _0, \end{aligned}$$
(62)

where after the derivation, all (external) fields are put to zero, while also the sources should be attributed their physical value (42), both manipulations denoted with \(\left. \right| _0\). \(\Gamma \) is the quantum effective action, which is subject to the same constrains (Ward identities) as its classical version \(\Sigma \). In particular, acting with \(\frac{\delta }{\delta \overline{\omega }_j^c}\) on

$$\begin{aligned}&\frac{\delta \Gamma }{\delta \omega _i^a(x)}+\partial _\mu ^x \frac{\delta \Gamma }{\delta N_\mu ^{ai}(x)}-gf^{abc}\overline{\omega }_i^b(x)\frac{\delta \Gamma }{\delta b^c(x)}\nonumber \\&\quad = -\hat{B} \partial _\mu ^x U_\mu ^{ai}(x)+D_\mu ^{ab,x}U_\mu ^{ib}(x) \end{aligned}$$
(63)

gives

$$\begin{aligned} \left. \frac{\delta ^2 \Gamma }{\delta \omega _i^a(x) \delta \overline{\omega }_j^c(y)}\right| _0=\left. -\partial _\mu ^x \frac{\delta ^2 \Gamma }{\delta N_\mu ^{ai}(x)\delta \overline{\omega }_j^c(y)}\right| _0, \end{aligned}$$
(64)

since we also have \(\frac{\delta \Gamma }{\delta b^a}=\partial _\mu A_\mu ^a\). The previous expression (64) can be further simplified, using the quantum version of the EOM (51):

$$\begin{aligned} \left. \frac{\delta ^2 \Gamma }{\delta \omega _i^a(x) \delta \overline{\omega }_j^c(y)}\right| _0&= -\left. \partial _\mu ^x \partial _\nu ^y\frac{\delta ^2 \Gamma }{\delta N_\mu ^{ai}(x)\delta U_\nu ^{cj}(y)}\right| _0\nonumber \\&+(1-\hat{B})\partial _\mu ^x \partial _\mu ^y\delta _{ij}\delta _{ac}\delta (x-y), \end{aligned}$$
(65)

while we keep the relations (42) in mind, together with global color invariance. Notice that, since there are no \((\omega ,\overline{c})\)-vertices, the contribution of the term \(\propto \frac{\delta ^2\Gamma }{\delta K \delta N}\) will vanish, as it corresponds to insertions of operators \(\sim c\overline{\omega }\) and no diagrams can be formed of this type.

Given \(\tilde{\Sigma }\) as in (37), (65) actually corresponds to

$$\begin{aligned} \left. \frac{\delta ^2 \Gamma }{\delta \omega _i^a(x) \delta \overline{\omega }_j^c(y)}\right| _0&= -\partial _\mu ^x \partial _\nu ^y \langle D_\mu ^{ak}\overline{\omega }_i^k(x) D_\nu ^{cm}\omega _j^m(y) \rangle _{1\mathrm{PI}} \nonumber \\&+\partial _\mu ^x \partial _\mu ^y\delta _{ij}\delta _{ac}\delta (x-y). \end{aligned}$$
(66)

We also need to reformulate the gap equation (34) a bit. Reconsider the relevant (physical) action (35). First we notice that the term \(\mathcal {T}\equiv - g f^{abc} \partial _\mu \overline{\omega }_i^a D_\mu ^{bd} c^d \varphi _i^c\) actually plays no role when computing quantities with Grassmann charge zero; indeed, there are no compensating \((\omega ,\overline{c})\)-vertices that would allow for the vertices of \(\mathcal {T}\) to generate any diagram. So, for practical purposes we can actually set \(\mathcal {T}=0\). Doing so, the action (35) enjoys the \(\varphi \leftrightarrow \overline{\varphi }\) symmetry; the extra \(\overline{\varphi }\varphi \partial A\)-term which comes about after interchanging the role of \(\varphi \) and \(\overline{\varphi }\) can be easily overcome by a shift in the multiplier \(b\)-field, or equivalently said, it vanishes on-shell because of the Landau gauge \(\partial A=0\). This observation can be used to reformulate the boundary condition (34)

$$\begin{aligned}&\left\langle gf^{abc}A_{\mu }^{a}\varphi _{\mu }^{bc}\right\rangle =\left\langle gf^{abc}A_{\mu }^{a}\overline{\varphi }_{\mu }^{bc}\right\rangle \nonumber \\&\quad =\frac{1}{2}\left\langle gf^{abc}A_{\mu }^{a}\left( \varphi _{\mu }^{bc}+\overline{\varphi }_{\mu }^{bc}\right) \right\rangle = d\hat{B}(N^2-1)\gamma ^2. \end{aligned}$$
(67)

Next, we notice that

$$\begin{aligned} d\hat{B}(N^2-1)\gamma ^2=\left\langle gf^{abc}A_{\mu }^{a}\overline{\varphi }_{\mu }^{bc}\right\rangle = \left\langle D_\mu ^{bc}\overline{\varphi }_\mu ^{bc}\right\rangle . \end{aligned}$$
(68)

Nextly, we consider the quantum Slavnov–Taylor identity,

$$\begin{aligned}&\int \mathrm {d}^4x\left( \frac{\delta \Gamma }{\delta K_\mu ^a}\frac{\delta \Gamma }{\delta A_\mu ^a}+\frac{\delta \Gamma }{\delta L^a}\frac{\delta \Gamma }{\delta c^a}+b^a\frac{\delta \Gamma }{\delta \overline{c}^a}+\overline{\varphi }_i^a\frac{\delta \Gamma }{\delta \overline{\omega }_i^a}\right. \nonumber \\&\quad \left. +\omega _i^a\frac{\delta \Gamma }{\delta \varphi _i^a}+M_\mu ^{ai}\frac{\delta \Gamma }{\delta U_\mu ^{ai}}+N_\mu ^{ai}\frac{\delta \Gamma }{\delta V_\mu ^{ai}}+R_\mu ^{ai}\frac{\delta \Gamma }{\delta T_\mu ^{ai}}\right) = 0\nonumber \\ \end{aligned}$$
(69)

and we act on it with \(\frac{\delta }{\delta N_\mu ^{ai}}\), then we set all external fields and sources (except for \(U\), \(M\), \(V\), \(N\), \(R\), \(T\)) to zero, thereby finding the identity

$$\begin{aligned}&\int \mathrm {d}^4 y \left( M_\nu ^{cj}(y)\frac{\delta ^2 \Gamma }{\delta N_\mu ^{ai}(x)\delta U_\nu ^{cj}(y)}\right) +\frac{\delta \Gamma }{\delta V_\mu ^{ai}(x)}\nonumber \\&\quad - \int \mathrm {d}^4 y\left( N_\nu ^{cj}(y)\frac{\delta ^2\Gamma }{\delta N_\mu ^{ai}(x)\delta V_\nu ^{cj}(y)}\right. \nonumber \\&\quad \left. +R_\nu ^{cj}(y)\frac{\delta ^2 \Gamma }{\delta N_\mu ^{ai}(x)T_\nu ^{cj}(y)}\right) = 0. \end{aligned}$$
(70)

Subsequently, we assign to the remaining sources their physical values. The contributions from the second line vanish, since \(N=0\) and because we can also drop the term \(\sim \) \(\frac{\delta ^2 \Gamma }{\delta N \delta T}\) since it corresponds to a correlation function of the type \(\langle \ldots \overline{\omega }c \rangle ^{1\mathrm{PI}}\), which is yet again zero since no diagrams can be formed because of the absence of \((\omega ,\overline{c})\) vertices. We also have

$$\begin{aligned} -\frac{\delta \Gamma }{\delta V_\mu ^{ai}}&= \left\langle D_\mu ^{ab}\overline{\varphi }_i^b\right\rangle ^{1\mathrm{PI}}+BM_\mu ^{ai}. \end{aligned}$$
(71)

Substituting (42) and (71) into (70), we obtain, taking into account the expression for \(\Sigma \) and the gap equation (68), while also keeping in mind the underlying global color and \(Q_i\)-charge invariance,

$$\begin{aligned} \int \mathrm {d}^d y \left\langle \left( D_\mu ^{ab}\overline{\omega }_\mu ^{ab}\right) _x \left( D_\nu ^{cd}\omega _\nu ^{cd}\right) _y\right\rangle ^{1\mathrm{PI}}=d\hat{B}(N^2-1). \end{aligned}$$
(72)

The \(1\mathrm{PI}\)-ness is of course coming from the fact that we worked with the \(1\mathrm{PI}\) generating functional \(\Gamma \). The limit \(d\rightarrow 4\) is implicitly understood.

In momentum space, Eq. (72) can be reformulated as a condition on the (\(1\mathrm{PI}\)) zero-momentum correlation function:

$$\begin{aligned} \lim _{p^2\rightarrow 0}\left\langle \left( D_\mu ^{ab}\overline{\omega }_\mu ^{ab}\right) \left( D_\nu ^{cd}\omega _\nu ^{cd}\right) _p\right\rangle ^{1\mathrm{PI}}=d\hat{B}(N^2-1). \end{aligned}$$
(73)

Specifically, parameterizing the \(1\mathrm{PI}\) correlation function as

$$\begin{aligned}&\left\langle \left( D_\mu ^{ab}\overline{\omega }_i^{b}\right) \left( D_\nu ^{cd}\omega _j^{d}\right) \right\rangle ^{1\mathrm{PI}}_p \nonumber \\&\quad =\delta _{ij}\delta _{ac}\left( \mathcal {F}_1(p^2)\delta _{\mu \nu }+\mathcal {F}_2(p^2)p_\mu p_\nu \right) , \end{aligned}$$
(74)

we get from (73)

$$\begin{aligned} \mathcal {F}_1(0)=\hat{B}, \end{aligned}$$
(75)

at least when assuming that \(\mathcal {F}_2(p^2)\) does not develop a dynamical pole at zero momentum. That the latter would be rather unlikely has been argued for by a Dyson–Schwinger analysis in [93].

Finally, inserting (66) into (61 76), we get

$$\begin{aligned} g(0)= \frac{1}{1-\mathcal {F}_1(0)}=\frac{1}{1-\hat{B}}=\frac{1}{B}, \end{aligned}$$
(76)

that is, the envisaged boundary condition.

4.4 Explicit one-loop analysis

It is instructive to check if the the boundary condition for the ghost propagator is indeed satisfied at the first order in the loop expansion. For this, we will need the vacuum energy and its behavior with respect to \(\gamma ^{2}\) and \(\hat{B}\). With this aim in mind, we only need the quadratic terms of (35) in order to obtain the vacuum energy and the ghost propagator. As some renormalization scheme will be used, the distinction between the bare and renormalized quantities will be explicitly written down. As usual the bare quantities will be identified by a subscript zero: as \(\phi _{0}\) for denoting fields and \(\lambda _{0}\) for denoting parameters. There are no specific labels for renormalized fields and parameters. The partition function accounting only for the quadratic terms of the proposed action can be recast in

$$\begin{aligned} Z_{\mathrm{quad}}&= \int [\mathcal{D}A[[\mathcal{D}c][\mathcal{D}\bar{c}] \exp \nonumber \\&\times \left[ -\int \frac{\mathrm {d}^{d}p}{(2\pi )^{d}}\, \left\{ \frac{1}{2}A^{a}_{\mu }(p) Q^{ab}_{\mu \nu }A^{b}_{\nu }(-p) \right. \right. \nonumber \\&\left. \left. - \bar{c}^{a}(p) p^{2} c^{b}(-p) - \gamma ^{4}\hat{B}d(N^{2}-1) \right\} \right] , \end{aligned}$$
(77)

where

$$\begin{aligned} Q^{ab}_{\mu \nu }&= \frac{1}{2}p^{2}\left[ \left( 1+\frac{2\gamma ^{4}g^{2}N}{p^{4}}\right) \delta _{\mu \nu }\right. \nonumber \\&\left. \quad -\left( 1-\frac{1}{\alpha }\right) \frac{p_{\mu }p_{\nu }}{p^{2}}\right] \delta ^{ab}. \end{aligned}$$
(78)

When the limit \(\alpha \rightarrow 0\) is taken at the end of a computation, the Landau gauge is recovered.

Therefore, as the (bare) vacuum energy is defined by

$$\begin{aligned} \mathrm{e}^{-VE} = Z, \end{aligned}$$
(79)

in the general case, we are able to find

$$\begin{aligned}&\mathrm{e}^{-VE} \nonumber \\&\quad = \exp \left\{ \gamma _{0}^{4}V\hat{B}_{0}d(N^{2}-1) - \frac{1}{2} \text {Tr} \ln Q^{ab}_{\mu \nu } + \text {Tr} \ln p^{2}\right\} ,\nonumber \\ \end{aligned}$$
(80)

which gives

$$\begin{aligned} E&= - \gamma _{0}^{4}\hat{B}_{0}d(N^{2}-1) +\frac{1}{2}(N^{2}-1)(d-1)\nonumber \\&\times \int \frac{\mathrm {d}^{d}p}{(2\pi )^{d}}\,\ln \left( \frac{p^{4}+2\gamma ^{4}Ng^{2}}{p^{2}}\right) -\int \frac{\mathrm {d}^{d}p}{(2\pi )^{d}}\,\ln p^{2}.\nonumber \\ \end{aligned}$$
(81)

For the ghost propagator, see Fig. 2, we need just to calculate the ghost form factor \(\langle \sigma (0)\rangle \), as can be seen form (21)–(23), keeping in mind that

Fig. 2
figure 2

Feynman diagram of ghost propagator at one-loop order

$$\begin{aligned} \left\langle A^{a}_{\mu }(p)A^{b}_{\nu }(-p)\right\rangle = \frac{p^{2}}{p^{4}+2\gamma ^{4}g^{2}N}P_{\mu \nu }\delta ^{ab}. \end{aligned}$$
(82)

Then we have

$$\begin{aligned} \langle \sigma (0)\rangle&= \frac{g_{0}^{2}N}{V(N^{2}-1)}\frac{1}{p^{2}} \int \frac{\mathrm {d}^{d}k}{(2\pi )^{d}}\, \frac{(p-k)_{\mu }p_{\nu }}{(p-k)^{2}}\frac{k^{2}}{k^{4}-2g^{2}N\gamma ^{4}}\nonumber \\&\times \left( \delta _{\mu \nu }-\frac{k_{\mu }k_{\nu }}{k^{2}}\right) , \end{aligned}$$
(83)

which can be reduced to

$$\begin{aligned} \langle \sigma (0)\rangle = \frac{g_{0}^{2}N}{V(N^{2}-1)}\frac{(d-1)}{d} \int \frac{\mathrm {d}^{d}k}{(2\pi )^{d}}\, \frac{1}{k^{4}-2g^{2}N\gamma ^{4}}. \end{aligned}$$
(84)

Making use of \(\overline{\text{ MS }}\) renormalization scheme after dimensional regularization, with \(d = 4-\varepsilon \) (\(\varepsilon \rightarrow 0\)) and with the result (48), we are led to the following renormalized quantities:

$$\begin{aligned} E=-\frac{2(N^2-1)}{g^2N}\hat{B}\lambda ^4+\frac{3(N^2-1)}{64\pi ^2}\lambda ^4\left( \frac{8}{3}-2\ln \frac{\lambda ^2}{\overline{\mu }^2}\right) , \end{aligned}$$
(85)

for the vacuum energy, and

$$\begin{aligned} \langle \sigma (0) \rangle = \frac{3}{4}\frac{g^{2}N}{(N^{2}-1)}\frac{1}{(4\pi )^{2}}\left( \frac{5}{6}-\ln \frac{\lambda ^{2}}{\bar{\mu }^{2}}\right) \end{aligned}$$
(86)

for the ghost form factor, upon setting \(\lambda ^4=2g^2N\beta \).

In order to evaluate (86) we need to fix the value of \(\lambda ^{2}\) (\(\sim \) \(\sqrt{\beta }\)), which is done by solving the gap equation,

$$\begin{aligned} 0&= \left. \frac{\partial E}{\partial \lambda ^2}\right| _{\lambda ^2=\lambda ^{*2}}= -\frac{4(N^2-1)}{g^2N}\hat{B}\lambda ^{*2}\nonumber \\&+\frac{3(N^2-1)}{32\pi ^2}\lambda ^{*2}\left( \frac{8}{3}\!-\!2\ln \frac{\lambda ^{*2}}{\overline{\mu }^2}\right) \!-\!\frac{3(N^2-1)}{32\pi ^2}\lambda ^{*2}=0,\nonumber \\ \end{aligned}$$
(87)

leading to the following vacuum energy and ghost form factor, respectively,

$$\begin{aligned} E=\frac{3(N^2-1)}{64\pi ^2}\lambda ^{*4} \end{aligned}$$
(88)

and

$$\begin{aligned} \langle \sigma (0) \rangle = \hat{B}. \end{aligned}$$
(89)

Therefore, the boundary condition over the ghost propagator at first-loop order was explicitly obtained from the proposed action (35).

Interestingly enough, the second derivative at the critical value is negative,

$$\begin{aligned} \left. \frac{\partial ^2 E}{(\partial \lambda ^2)^2}\right| _{\lambda ^2=\lambda ^{*2}}=-\frac{3(N^2-1)}{16\pi ^2}, \end{aligned}$$
(90)

meaning that we are dealing with a local maximum rather than minimum, putting in a different perspective the comments made below (28). Notice that the general argumentation made there and the explicit one-loop verification are not necessarily at odds with each other. An important subtlety not taken into account below (28) is, besides that in practice we work in a formal series expansion, one must renormalize, which means infinities are subtracted. It is immediate to realize that a positive infinity can become negative after the renormalization procedure. The derivation of the action and gap equation must be understood at the formal and bare level, in which case there is no problem with negative variance \(\langle (H-\langle H \rangle )^2 \rangle \).

The gap equation (87) can be explicitly solved for by summing the potentially large logarithms. This is achieved by setting \(\overline{\mu }^2=\lambda ^{*2}\) (corresponding to use the renormalization group invariance of the effective action), giving

(91)

where the one-loop coefficient of the \(\beta \)-function reads \(\beta _0=\frac{11}{3}\frac{N}{16\pi ^2}\). Henceforth,

(92)

Having found the expression for the vacuum energy in terms of the parameter \(\hat{B}\), one might wonder if it would be possible to find a specific value for \(\hat{B}\) by using it as a variational parameter. In the CJT formalism [94], the two-point functions are to be determined by minimizing the effective action w.r.t. those.Footnote 12 However, this frequently amounts to a tremendously difficult task, in which case it can be fortuitous to model the propagators in terms of a (few) parameter(s), leading to algebraic equations for these parameters from the minimization of the effective action. From this perspective, we can indeed employ \(\hat{B}\) as a variational parameter associated to a specific ghost propagator and discuss whether some \(\hat{B}\) minimizes the corresponding vacuum functional.

In a one-loop approximation, given the result (92), it is immediately clear that \(\frac{\partial E}{\partial \hat{B}}=0\) has no sensible solutions in terms of \(\hat{B}\). This is not a coincidence, as we can provide an exact argument valid at any order. Writing the vacuum functional as

$$\begin{aligned} E(\lambda ^2)=\tilde{E}(\lambda ^2)-\frac{2(N^2-1)}{g^2N}\hat{B}\lambda ^4, \end{aligned}$$
(93)

we need to solve, for each value of \(\hat{B}\),

$$\begin{aligned} 0=\left. \frac{\partial E}{\partial \lambda ^2}\right| _{\lambda ^2=\lambda ^{*2}}= -\frac{4(N^2-1)}{g^2N}\hat{B} \lambda ^{*2}+\left. \frac{\partial \tilde{E}_{vac}}{\partial \lambda ^2}\right| _{\lambda ^2=\lambda ^{*2}}.\nonumber \\ \end{aligned}$$
(94)

The vacuum energy, \(E(\lambda ^{*2})\), must then depend minimally on \(\hat{B}\). The condition can be written as

$$\begin{aligned} 0=\left. \frac{\partial E}{\partial \hat{B}}\right| _{\lambda ^2=\lambda ^{*2}}&= \left. \frac{\partial \tilde{E}}{\partial \lambda ^{2}}\right| _{\lambda ^2=\lambda ^{*2}}\frac{\partial \lambda ^{*2}}{\partial \hat{B}}-\frac{4(N^2-1)}{g^2N}\nonumber \\&\times \hat{B}\lambda ^{*2}\frac{\partial \lambda ^{*2}}{\partial \hat{B}}-\frac{2(N^2-1)}{g^2N}\lambda ^{*4}. \end{aligned}$$
(95)

Upon using the supplementary condition (94), the expression (95) reduces to

$$\begin{aligned} -\frac{2(N^2-1)}{g^2N}\lambda ^{*4}=0, \end{aligned}$$
(96)

requiring that \(\lambda ^{*2}=0\), which is not allowed, otherwise the gauge field region of integration will not be restricted to the first Gribov region and, thus, the ambiguity problem in gauge fixing will be present. This means that using our formalism, there is no optimal value for \(\hat{B}\) to be found that would minimize the vacuum energy, at least when restricting the analysis to a specific order in perturbation theory.

5 Discussion

In the first part of this paper, we have tried to take into account at the level of the path integration the recently derived Cucchieri–Mendes bounds on the ghost propagator in the presence of an external gauge field that belongs to the Gribov region \(\Omega \). Our result indicates that this bound might not leave its footprints in the continuum restriction of the gauge field path integral, since the integral gets its main contribution from gauge fields close to the boundary where the (stricter) bound of Cucchieri–Mendes is saturated. Our result is based on an infinite-volume continuum formulation, so we cannot say much about what happens exactly when a genuine discrete setting (as on a lattice, the setting of [43]) would be used, e.g. how “fast” the boundary of the Gribov region is approached or how this would depend on a variable lattice size/volume.

In the second part of the paper, we have provided a positive answer to the query put forward by Maas [45, 46] by constructing an action principle that allows one to give a prescribed value to the zero-momentum ghost form factor and this without jeopardizing the renormalizability. We have also given a formal proof, based on the Ward identities of the theory, that the boundary condition is met. One issue we did not touch upon is that we did not show that our implementation of the restriction corresponds to a stable theory. Likely, following [811], the vacuum will be unstable against the condensation of \(d=2\) operators which will drastically effect the properties of e.g. the propagators, given the underlying vacuum will get nonperturbatively changed. If this happens, the Ward identities of the new theory will be different, and the proof that the boundary condition is realized, will no longer apply. Said otherwise, the ghost form factor at zero momentum, \(g(0)\), will attain a value that will be not predetermined, but needs to be computed in one or another approximation. This scenario would correspond to a quantum modification of the fact that the path integral is geometrically dominated by configurations near the boundary \(\partial \Omega \) of the Gribov region.

An important ingredient in the refinement of Gribov–Zwanziger like theories is the condensation of the \(d=2\) operator \(\overline{\varphi }_\mu ^{ab}\varphi _\mu ^{ab}-\overline{\omega }_\mu ^{ab}\omega _\mu ^{ab}=-s(\overline{\omega }_\mu ^{ab}\varphi _\mu ^{ab})\). This operator can be coupled to the action \(\Sigma \) (43), via the addition of \(-\int \mathrm {d}^dx J(\overline{\varphi }_\mu ^{ab}\varphi _\mu ^{ab}-\overline{\omega }_\mu ^{ab}\omega _\mu ^{ab})\) with \(J\) a local BRST invariant source, \(sJ=0\). The Ward identities used will thus at most get modified by linear breakings in \(J\), so they remain useful at the quantum level. The main obstacle, however, is that we lose the connection (60) between the Faddeev–Popov ghost propagator and the auxiliary Gribov–Zwanziger ghost propagator, since the latter acquires a dynamical mass due to \(J\) (more precisely, due to the ensuing condensation; see [810]). So, even if we would be able to prove something about the \(\langle \overline{\omega }\omega \rangle \) propagator at zero momentum, we do not longer have a connection to the propagator of interest, \(\langle \overline{c} c \rangle \). Also in functional formalisms, a similar caveat might interfere. One can write down the quantum (Dyson–Schwinger) equations of motion based on the action (33). Since the Ward identities should in principle be obeyed by the nonperturbative \(n\)-point functions,Footnote 13 the boundary condition of the ghost propagator should remain valid. Though, if composite operators, as \(\overline{\varphi }_\mu ^{ab}\varphi _\mu ^{ab}-\overline{\omega }_\mu ^{ab}\omega _\mu ^{ab}\) can condense, their dynamics can be included in the functional formalism as well, e.g. by applying the CJT method [94]. Similarly, the Ward identities of the CJT action will be adapted, changing the conclusion about the ghost boundary condition.

Lastly, an important feature of our formalism, completely analogous to the standard Gribov–Zwanziger formalism, is that it comes with a soft breaking of the BRST. Indeed, taking for instance the \(s\)-variation of the action (33), we get

$$\begin{aligned} sS&= \int d^4x\left( \sqrt{\beta ^*} gf^{abc}D_\mu ^{am}c^m\left( \varphi _\mu ^{bc}+\overline{\varphi }_\mu ^{bc}\right) \right. \nonumber \\&-\left. \sqrt{\beta ^*}gf^{abc}A_\mu ^a\omega _\mu ^{bc}\right) \end{aligned}$$
(97)

which is clearly nonzero if \(\beta ^*\ne 0\). In retrospect, this also indicates we cannot give \(\hat{B}\) the same meaning as a genuine gauge parameter. Clearly, physical quantities will depend on it. A first explicit example is the vacuum energy that varies with \(\hat{B}\). If we define physical operators as the quantum extensions of the classical gauge invariant operators (see e.g. [95]), their expectation values will in general depend on \(\hat{B}\). First of all, \(\hat{B}\) is not coupled to a BRST-exact expressionFootnote 14 and even if it would be, the BRST symmetry is broken, so the standard argument that the expectation value of BRST-exact operators nullifies would anyhow no longer apply.Footnote 15

As the BRST breaking is proportional to a lower dimensional operator, it is a soft breaking and hence controllable at the UV level. This is at the core of the renormalization properties of the GZ-like theories. The precise physical consequences of such a BRST breaking are not fully understood yet [96], despite a series of various recent efforts [97102]. Here, we only wish to point out that in a recent numerical effort, further lattice backup has been provided in favor of the BRST breaking by a restriction to the first Gribov region [103]. More work is under way concerning the search for a (lattice verifiable) signal of the BRST breaking [104, 105].