Landau gauge Yang-Mills propagators in the complex momentum plane

We calculate the dressed gluon and ghost propagators of Landau gauge Yang-Mills theory in the complex momentum plane from their Dyson-Schwinger equations. To this end, we develop techniques for a direct calculation such that no mathematically ill-posed inverse problem needs to be solved. We provide a detailed account of the employed ray technique and discuss a range of tools to monitor the stability of the numerical calculation. Within a truncation employing model ansaetze for the three-point vertices and neglecting effects due to four-point functions, we find a singularity in the gluon propagator in the second quadrant of the complex $p^2$-plane. Although the location of this singularity turns out to be strongly dependent on the model for the three-gluon vertex, it occurs always at complex momenta for the range of models considered.


I. INTRODUCTION
There are at least two reasons why the analytic structure of Yang-Mills propagators, viz., of the ghost and the gluon propagators, are of great interest. First, there are direct connections to fundamental properties and problems of the theory such as color confinement, the associated construction of an asymptotic state space in terms of gauge invariant and colorless states, and the question whether BRST-symmetry is broken non-perturbatively or not. Second, on a practical level, the analytic structure of the propagators plays an important role in the calculation of all properties of bound states, not least glueballs, in the functional framework of Dyson-Schwinger and Bethe-Salpeter equations. It is also natural to assume that these two issues are related to each other. Consequently, the topic received a lot of attention over the years. In the past century, based on studies of the gauge fixing problem, Gribov [1] and Zwanziger [2] suggested an explicit expression for the gluon propagator with complex conjugate poles at purely imaginary squared Euclidean momenta p 2 . Stingl [3] provided a generalization of this ansatz by shifting the complex conjugate poles to general complex momentum squares in the negative half-plane. The refined Gribov-Zwanziger framework formulated later leads also to conjugate poles located at complex momenta [4]. An alternative form with a branch cut structure for real and timelike squared momenta was proposed in Ref. [5]. In recent years, research focused in addition on general properties of the spectral function of the gluon, which in turn restricts its analytic structure [6][7][8][9][10]. There are in principle two different strategies to extract the analytical properties of the gluon from explicit re-sults of non-perturbative approaches such as lattice or functional methods. Lattice Yang-Mills theory generically delivers results for positive and real (i.e. spacelike) momenta. Thus, reconstruction methods have to be employed to study the analytic continuation into the complex momentum plane. In this respect, many of the above mentioned explicit forms have been used as trial functions to describe lattice data at real and spacelike p 2 [11][12][13][14]. In addition, reconstruction algorithms like the Bayesian spectral reconstruction method, the Tikonov regularization or the Schlessinger point method have been used [15][16][17]. Of course, these reconstruction methods can be applied equally well to solutions from functional methods, i.e. either Dyson-Schwinger equations or the functional renormalization group [8,17,18]. In addition, such functions can also be used to analytically continue results (instead of correlation functions) obtained from Euclidean input to the physical momentum regime. This was successfully realized for the calculation of (pseudo)scalar glueballs [19], where the availability of Euclidean input from a self-contained calculation [20] led to results in quantitative agreement with lattice results [21][22][23]. One of the advantages of the functional approach, however, is that direct calculations at timelike momenta are possible. This property is exploited routinely in the calculation of spectra and properties of bound states, see e.g. [24,25] for reviews. For the gluon propagator, a first explicit calculation was discussed in [26] using a particular technique (the 'ray method') that has been developed in the context of QED in three dimensions [27]. It was also used for several other purposes since then, e.g., [5,[28][29][30][31][32]. Subsequently, other techniques for the gluon propagator, in particular a direct solution on a momentum grid in the complex plane, were also explored [33]. In this work, we expand upon and refine previous work using the ray technique [26]. We improve the numerical stability of the method and introduce a number of tools to monitor the reliability of the obtained results on a step by step basis when probing the complex momentum plane further towards the timelike region. As a result, we are able to resolve an analytical structure in the second quadrant which was not seen in [26]. We update and discuss the corresponding results for the gluon and the ghost propagators. The remainder of this article is organized as follows: In Sec. II, we explain the underlying idea of the ray technique. The setup of our calculations is discussed in Sec. III, and the results are presented in Sec. IV. In this section, also several tests are introduced and applied to the results. We close with a summary in Sec. V. Computational details, the reconstruction from arbitrary rays and the employed three-gluon vertex models are explained in appendices.

II. THE RAY TECHNIQUE
The exact set of coupled Dyson-Schwinger equations (DSEs) for the gluon and ghost propagators in Landau gauge Yang-Mills theory are displayed in Fig. 1. In Landau gauge, the ghost and gluon propagators, D G and D µν , are given by where color factors have been suppressed. Their DSEs feature nonperturbative one-and two-loop diagrams on the right-hand side, which all share one property: If the momentum variable p 2 that enters the diagrams from the outside is complex, poles and branch cuts in the various integrands appear. In principle, this is the reason why a direct solution of these equations on a complex momentum grid is extremely dangerous if not prohibitive, since it automatically implies integration across cuts. In the quark sector of quantum chromodynamics, the situation is somewhat alleviated by the quark mass, which modifies the location of these cuts [29,34] and allows the calculation in a restricted momentum region. In practical calcu-lations using rainbow-ladder type models, it depends on the type of the model whether cuts are absent [35,36], small [37] or potentially relevant on a quantitative basis [38]. For the gluonic system, such a rainbow-like truncation was employed in Ref. [33]. However, due to the structure of the integrals in the gluon propagator DSE, this breaks the self-consistency of the equations, because the propagators one would like to solve for are also contained implicitly in the models. If we want to maintain self-consistency, the appearance and proper treatment of cuts in the integrands seem unavoidable [26]. We therefore need a different strategy [27], which we call the 'ray technique'. We illustrate the basics of the ray technique using a simple massless scalar model with a cubic interaction 1 . The corresponding selfenergy diagram in the DSE for the scalar propagator has the same structure as those we consider in Yang-Mills theory but without the complications of Lorentz tensors. The massless propagator is described by with x = p 2 and Z(x) its dressing function. The perturbative one-loop self-energy is given by where y = q 2 , z = (p + q) 2 = x + y + 2 √ x y cos θ and constant factors were dropped in the second line for brevity. This notation is kept throughout this paper, i.e., the external momentum squared is x, the internal one squared y and the squared combined momentum z. We regularize the integral by the O(d) symmetric UV cutoff Λ for the radial part. For now, we also keep the dimension d general.
Clearly, the integrand is singular for z = 0. For fixed external momentum x, this singularity is located at y = x exp(±2 i θ). The angular integration over θ then leads to a branch cut in the y-plane with the endpoints at y = x. As is visualized in Fig. 2, the branch cut lies on a circle with |y| = |x|. There is only one point at arg x = arg y where the cut is open. For any non-real or negative x, the usual integration path of y along the positive real axis is now forbidden, since it would cross the cut. To avoid this problem, one needs to deform the integration path from y = 0 to y = Λ 2 such that it goes through the opening at y = x in between. However, it depends on the dimension (and in the most general case on the details of the integral kernels and the dressing functions) whether the opening is suitably finite and the path is save. The crucial factor here is the behavior of the integrand at the boundaries of the angular integration. In our example, the perturbative treatment of the scalar theory, singularities for |y| = |x| appear in d = 2 dimensions which is just a manifestation of perturbation theory being ill-defined for the massless theory in two dimensions. In higher dimensions, the path deformation is possible. For Yang-Mills theory, the situation is more complicated 2 . Let us first have a closer look at the DSE for the ghost propagator which contains the ghost-gluon vertex as only quantity beyond the propagators. It is a very well known object with only one dressing function D Acc,T contributing to the integrand due to the transversality of the gluon propagator in Landau gauge. The behavior of the integrand of the angular integral is then controlled by a momentum dependent kinematic kernel times the dressing functions of the vertex and, depending on the momentum routing, the dressing of the ghost or the gluon propagator. For the latter situation we obtain where the extra factor (sin θ) 2 in the kernel stems from the contraction of the gluon with the two vertices. This angular integral is finite at |y| = |x| in any dimension provided the vertex dressing function D Acc,T does not develop a strong singularity for z = 0 that overcompensates for the kernel and the well-known infrared behavior Z(z) → 0 of the gluon dressing function. The information we have about the ghost-gluon vertex does not support the existence of such a problematic singularity, e.g., [20,[40][41][42][43][44][45][46][47][48][49][50][51][52]. A similar situation arises for the ghost loop in the gluon propagator DSE. For the gluon loop, however, terms proportional to 1/z 2 appear (see, e.g., [53] for explicit expressions). Fortunately, these are countered by the presence of at least one factor of Z(z) → 0 in the integrand. Thus, again, provided the three-gluon vertex does not develop a strong singularity at z = 0, the integral is finite and the path deformation works. It should be stressed here that the crucial momentum variable z affects only one leg of the three-gluon vertex. The relevant divergence structure is thus that of one momentum going to zero and not the global IR behavior of the three-gluon vertex. Luckily, the former divergences were found to be only weak [42][43][44].
One can show as well that the path deformation works for the two-loop diagrams. Since below we will deal with a truncation involving one-loop diagrams only, we refrain from going into detail here and refer to Ref. [52] for details on their structure. We furthermore wish to emphasize that the problem of potential singularities at |y| = |x| is absent for massive propagators such as quarks. Due to the finite mass one then encounters an opening of the branch cut of finite size [29,34].
Additional problems may be encountered if dressing functions present under the integral develop poles or branch cuts at momenta other than z = 0 probed by the integration. In principle, this leads to additional constraints on the integration contour. A typical case is a pole in a propagator D(x) at complex momentum x 0 . If this propagator only appears under the radial integral, then the path may either be deformed around the pole, or a corresponding residue needs to be taken into account. The situation is worse if the propagator also depends on the angle, i.e. if we have D(z) or more general momentum arguments. The singularity condition is then z = x 0 and corresponds to the case of a massive propagator with mass m = −x 0 . Thus in principle it can be dealt with analogously, as already discussed above and in [29,34]. In practice, however, this solution is very hard to implement in a self-consistent way in numerically demanding situations such as the coupled system of ghost and gluon propagators.
Setting these potential additional problems aside for the moment, it remains to be discussed how the integration contour is chosen in practice. A typical integration path is shown in Fig. 3. The integration contour runs along a radial 'ray' and is then continued to the cutoff Λ 2 by a second curve which we call 'arc'. The precise form of the arc is not so relevant in practice. All details concerning the numerical implementation of the ray technique are discussed in Appendix A.

III. TRUNCATION AND RENORMALIZATION OF THE GLUON AND GHOST PROPAGATOR DSES
In this work, we are primarily interested in a conceptual study of the gluon propagator at complex momenta. Thus, although the two-loop diagrams in Fig. 1 are quantitatively important on a 20 % level [20,54], we neglect them to avoid substantial technical complications. The tadpole diagram is dropped as well because it vanishes in the renormalization we employ. This truncation leaves us with the system depicted in Fig. 4. This system is closed once the dressed ghost-gluon and three-gluon vertices are known. It is well known [20,[40][41][42][43][44][45][46][47][48][49][50][51][52]55] that the ghost-gluon vertex only receives small nonperturbative corrections and therefor in our conceptual study we take it as bare. For the three-gluon vertex we tested several models that will be described in the next subsection. The second part of this section describes the renormalization procedure.

A. The three-gluon vertex
The three-gluon vertex was studied in various approaches, ranging from lattice simulations [56][57][58][59][60][61] to effective models using a massive gluon propagator [47] to functional equations [20,42,43,49,[62][63][64][65][66]. From these studies, the nonperturbative properties for spacelike momenta are quite well understood. In particular, a suppression at intermediate momenta is seen, and the dressing of the tree-level tensor structure even becomes negative. However, the zero crossing is at rather low momenta which are difficult to reach in lattice calculations.
Here we use a model for the three-gluon vertex which also includes a term that restores the correct renormalization group behavior of the gluon propagator. Note that this is only necessary because we discarded the two-loop terms in the gluon propagator DSE. For more elaborate truncations, the correct renormalization group behavior is obtained automatically, see Refs. [20,52,67]. The vertex model is restricted to the tree-level tensor and parameterized as where Γ AAA,(0) µνρ (p, q, r) is the Lorentz tensor of the treelevel vertex. We tested various forms for C AAA (p 2 , q 2 , r 2 ) which are detailed in Appendix C. As the results for the gluon and ghost propagators are qualitatively very similar for these models, we choose one representative for illustration in plots. This model reads Z 1 is the renormalization constant of the three-gluon vertex, δ = −9/44 the anomalous dimension of the ghost propagator and a a parameter that determines the IR behavior of the model. Note that in the UV a drops out due to the scaling relation 1+2δ+γ = 0 of the anomalous dimensions δ and γ of the ghost and gluon propagators, respectively. The vertex model is a reparametrization of the model introduced in Ref. [45] without IR part and corresponds to a Bose symmetrized version of the model from Ref. [53]. The non-Bose symmetric version is recovered by replacing the dressings as G(p 2 ) 2 → G(y)G(z) and Z(p 2 ) 2 → Z(y)Z(z). This version was used in Ref. [26] and also tested here. Again, though, we did not find any qualitative differences. The IR behavior of this model is determined by the parameter a. Originally, a = 3δ was used which makes the expression IR finite for the scaling solution. Here, we solve for a decoupling solution for which C AAA 1 (x, y, z) is IR divergent with a = 3δ. 3 One way to circumvent this is to use a = −1 instead. This removes the gluon dressing function from the model. An alternative way is to modify the momentum argument by adding a small scale Λ 2 IR . We do so by setting In the plots shown in Sec. IV, we used a = 3δ with Λ 2 IR = 0. However, we checked that the results from the two methods are qualitatively the same.
In this model, the analytic behavior of the three-gluon vertex is completely determined by the analytic structure of the propagator dressing functions. As we found that for the given truncation the gluon propagator seems to have a singular point in the complex plane, we wanted to remove its influence on the vertex. The corresponding models are discussed in Appendix C. While these adaptations do have quantitative effects, all qualitative aspects remain the same. Thus, we conclude that our qualitative results are robust against changes within a large class of three-gluon vertex models, but we stress that our analysis only holds for this class.

B. Renormalization
The propagator DSEs are renormalized via a momentum subtraction scheme. For the ghost propagator DSE this leads to G(s) denotes the ghost dressing function at an (infrared) subtraction scale s and Σ G (p 2 ) is the ghost selfenergy. The value chosen for G(s) selects a particular solution from a one-parameter family of possible ones, for details see [70]. To be able to perform the subtraction ray by ray we need to analytically continue the value of G(s) from ray to ray. We do so using the Cauchy-Riemann condition, see Eq. (10) below and Appendix A for details. For the gluon propagator DSE, the subtraction point needs to be chosen at large momenta for numerical stability. The hard UV cutoff employed in our calculations breaks gauge covariance and leads to additional quadratic divergences in the gluon propagator DSE. Several methods to remove them exist, see Ref. [75] and references therein for details. Such a procedure can entail ambiguities and should thus be considered as part of the model input. However, a unique subtraction might be possible with more elaborate truncations [20]. Since the present truncation scheme is constructed for qualitative tests, we use a simple but effective subtraction scheme that modifies the integrand of the gluon loop appropriately [53,76]. We also tested other methods. For one, instead of subtracting the overall divergence in the gluon loop, we split the subtraction between ghost and gluon loops [45]. For another one, we used a second renormalization condition [67,77] which already proved very useful elsewhere [20]. It turns that out the second renormalization condition can be chosen such that both procedures lead to almost similar results. Finally, it remains to set the physical scale of our results. We do so by matching the gluon dressing function to corresponding lattice results [78]. Thus, we inherit the scale setting used on the lattice.

A. Baseline setup
In this section we first present results for what we call 'baseline setup'. It is defined by fixed values for the renormalization conditions, the subtraction method of quadratic divergences in the gluon loop and the threegluon vertex model from Eq. (7) with a = 3δ. More details including computational parameters are explained in Appendix A. Variations and tests of this setup are presented in the subsequent subsections. The gluon dressing function is shown in Fig. 5 as a function of the radial and angular parts of the complex variable x = p 2 =p 2 e iθ . In the first quadrant, the calculation works without problems. However, in the second quadrant, the bump in the gluon dressing function atp 2 ≈ 0.5 starts to rise appreciably, however, without becoming singular. Beyond a certain value of θ, it flattens again and remains finite until the negative squared momentum axis. This behavior was seen before in a slightly different setup and with less precise numerics [26] and was interpreted as the gluon being regular at complex momenta. However, with the improved numerical treatment followed in this work, we clearly observe oscillations in the solution, which may hint towards numerical artifacts. These signals and the strong rise could also indicate that we may have hit a singular point beyond which the ray technique is no longer applicable. We tried to take this finding into account by modifying the integration path appropriately (cf. the discussion in section II), but then we loose the advantage of the ray technique that we do not need to know the dressing function in unknown regions, and we did not succeed in improving the results in this way. Thus, from this plot we concluded that our results probably cannot be trusted beyond that point and that we may have hit a singularity. We corroborated this conclusion further using various tests that will be detailed below in section IV B.
Our results for the ghost dressing function in the baseline setup are shown in Fig. 6. In contrast to the gluon propagator, we do not see any drastic changes in the dressing functions. The real part of the ghost is smooth throughout the complex momentum plane, whereas the imaginary part develops a negative bump at the same scale at which the gluon rises massively. However, since the ghost and the gluon equations are directly coupled, the results for the ghost dressing should only be considered trustworthy up to the location of the potentially singular point in the gluon dressing function. The reason for the appearance of the potential singularity is not clear. In particular, it cannot be understood in simple terms similar to the Cutkosky rules [79,80] which allow to determine the position of a branch cut from the masses of the propagators in a Feynman diagram. In the language of contour deformation, a branch cut in the external momentum arises when the integration path cannot be deformed continuously for two values of p 2 , see Ref. [34] for details. This is, however, not the case here.

B. Variations
Since the origin of the potential singularity thus seems to be dynamic, we tried to vary our setup to investigate the influence on the existence and position of the singular point. We attempted the following variations of the baseline setup: • We tried several modifications of the three-gluon vertex model. The corresponding expressions are listed in Appendix C. In addition, we varied the parameter a in the model between 8δ ≈ −1.61 and 3δ ≈ −0.61. The effect of this is discussed below.
• We tested alternative methods to subtract quadratic divergences as discussed in Sec. III B.
• We varied the renormalization condition G(0) of the ghost propagator to obtain different decoupling solutions.
Most variations did not lead to a qualitative change of our results and we do not discuss them further. Different values for G(0) only led to small quantitative changes in the position of the singularity. A larger effect was observed from the employed three-gluon vertex model. In particular, the position of the potentially singular point in the complex momentum plane for the gluon propagator depends on the details of the vertex model. Varying the parameter a, cf. Fig. 7, we observed that the point moves closer to the real axis (and at the same time closer to the origin) when the parameter a was lowered. Lowering a as far as a = 8δ ≈ −1.61, we were able to move the potential singularity almost onto the real axis. This, however, happens on the expense that gluon dressing function at real spacelike momenta becomes unrealistically flat. We therefore did not lower a further.

C. Tests
A propagator can be calculated in the complete complex plane from the spectral density σ(s) via if no poles at complex momenta exist, see also Appendix B. The inverse task of extracting the spectral density from the propagator given at Euclidean momenta is an ill-posed inverse problem [81]. This is reflected in the necessity of some form of bias in these methods and a large sensitivity of the results on the precision of the input. The direct calculation performed here, on the other hand, does not have these intrinsic problems. Rather, the main challenges are of numeric nature, see Appendix A.
In particular, the global nature of analyticity can be problematic, viz., the analytic properties of a function are encoded in the behavior of the function on any region of the complex plane. Hence, the propagation of errors has to be under control and it is important to check that the numeric calculation does not interfere with analyticity. In this section, we describe several possibilities for such checks and apply them to our results. One direct possibility for such a check is provided by Cauchy's integral formula. We use it to reconstruct the propagators for real and spacelike momenta from any ray and monitor the quality of the reconstruction. The details for this procedure are listed in Appendix B. We performed this reconstruction on all rays used for the calculation of the gluon propagator and ghost dressing functions. For the baseline setup, the reconstructed functions are shown in Fig. 8. For guidance, we also plot the Euclidean result (in red at θ = 0). The reconstruction works very well in the first quadrant (θ ≤ π/2) and some way into the second quadrant until it fails on the ray marked with a thick black line. Beyond this ray, it fails completely. Again, this indicates the appearance of a singularity in the complex plane at or around the ray in black. Another possibility for a check uses the Cauchy-Riemann equations which relate the real and imaginary parts of analytic functions. In polar coordinates they can be written as where x = r e i θ . We can use this to test the analyticity of our results by monitoring We approximate the derivatives by finite differences which leads to its numeric deviations on its own. The results are shown in Fig. 9. Again, we find small deviations between the calculated and the Cauchy-Riemann extrapolated results for ghost and gluon propagators (cf. the scale of the z-axis) up to the singular point in the second quadrant beyond which the Cauchy-Riemann equations are clearly no longer fulfilled. The two ridges are artifacts from splitting the grid for the dressing functions as discussed in Appendix A.
Having collected ample evidence that our results are analytic until we hit a potential singularity in the second quadrant, we turn the situation around and test the reliability of extrapolating the results from the spacelike axis using the analytic continuation from a Padé approximant. Specifically, we use the Schlessinger point method [82] which provides a pointwise exact description of the data. As input, we use a random subset of momentum points on the positive and real momentum axis. Note that this test has been performed already in Ref. [17] for one of the truncations that we also used here. There, indeed a singularity in the second quadrant was found. Here, we check whether this finding persists for all truncations considered in this work and, even more important, whether the results from the Schlessinger point method agree quantitatively with the ex- plicit results from solving the DSEs in the complex momentum plane. We compare the direct calculation and the Schlessinger point method by plotting the ratio where f SPM (x) is the function as obtained from the Schlessinger point method. This is shown in Fig. 10. The plots confirm that up to the potentially singular point the two methods agree quite well. Next, we compare the positions of the singular points found by the two methods for the gluon propagator. To this end, we follow a simple procedure. For a given subset of input momentum points on the real and positive momentum axis, we calculate the pole positions analytically from the coefficients of the Schlessinger point method. Poles with small residues are discarded as artifacts. For a final estimate of pole positions, we sample several subsets of the Euclidean data. This method is not as elaborate as the one introduced in Ref. [17], where the sample of input points is optimized based on the quality of the reconstruction from the spec- Thus, the combined evidence of all methods clearly points towards a singularity in the gluon dressing function at complex momenta. We observed the same quantitative agreement between the pole location predicted by the Schlessinger method and the explicit results from the ray method for all truncations studied in this work.

V. SUMMARY
In this work, we studied the analytic structure of the gluon and ghost propagators of Landau gauge Yang-Mills theory from a coupled set of Dyson-Schwinger equations.
Using the ray technique, we were able to solve the equations for a region of complex squared momenta that extended well into the second quadrant. We continuously checked our calculation by a variety of other methods, namely by (i) reconstruction algorithms using our solutions to reconstruct the propagators on the spacelike real momentum axis, (ii) the Cauchy-Riemann equations, and (iii) the Schlessinger point method which provides rational functions for the analytic continuation that can be compared with our explicit results. All of these methods agree very well up to a certain ray in the second quadrant. At this point, we encountered a steep rise in the gluon dressing function at the same location where the Schlessinger method predicts a pole. Thus, the combined evidence of both methods strongly suggests the presence of a non-analytic structure in the complex plane. Due to the less precise numerics available at the time, this structure was not recognized as a singularity in Ref. [26].
We studied the properties of this singularity and noted that within our truncation, the details of the model for the three-gluon vertex are most relevant for its location, whereas other technical details such as the renormalization procedure matter much less.
It will be important in the future to further check the dependence of the analytic structure on the gluon interaction. In this respect, it is highly relevant to improve this calculation with better input for the three-and fourgluon vertices, e.g. by using explicit input from solutions of their respective DSEs.

VI. ACKNOWLEDGMENTS
We are grateful for intense discussions with Arno Tripolt in the early stages of this work. We furthermore acknowledge fruitful interactions with Gernot Eichmann, Joannis Papavassiliou and Jan Pawlowski. This work was supported by the Helmholtz Research Academy Hesse for FAIR, by the DFG (German Research Foundation) grant FI 970/11-1, and by the BMBF under contract No. 05P18RGFP1.
Appendix A: Techniques used to solve the DSEs numerically To extract a high quality numerical solution of the coupled system of DSEs for the gluon dressing function Z(p 2 ) and the ghost dressing G(p 2 ) in the complex p 2 momentum plane, we employ a variety of numerical tools that are described in some detail in the following. Some of these have been already used in a previous publication [26].

The grid
As explained in section II, we solve the DSEs on a grid of 'rays' and 'arcs'. Each ray extends radially outwards from the origin to a fixed momentum cutoff at p 2 = e iθ Λ 2 1 , where θ denotes the angle between the ray and the positive real axis. From this point on, the 'arc' connects the ray with the real axis along a path given by with t ∈ [0, 1]. Thus, we have to deal with two different cutoff scales, Λ 1 and Λ, which are chosen to be close to each other, see Tab. I. For the grid of rays we typically use 180 rays that cover the complex plane from θ = 0 to π− . We also tested using 360 rays or 90 rays instead but did not find any significant influence on the final results. On the real axis, the corresponding ray and arc are of course collinear and merge into a straight line integration path.

Representation of the dressing functions
The integration paths in the DSEs are along the rays and along the arcs. Under the integrals of the DSEs, we encounter two different types of arguments in the dressing functions Z and G: On the one hand, there is the integration momentum denoted by q 2 = y. On the other hand, there is the squared difference between external momentum p and integration momentum q denoted by (q − p) 2 = z. The external momenta p 2 = x are distributed over the rays, while the integration momenta can be on the rays and the arcs. If y is on a ray, the squared differences, z, are also on the rays or their extensions. When y is on an arc, however, z can also take values elsewhere in the complex plane, see Fig. 11 for two examples.
To carry out the integrations on the right hand side of the DSE, we need the dressing functions Z and G on all these points. In the following we explain in detail how we manage this. Let us first deal with the rays. On each ray, we represent the real part and the imaginary part of Z and G separately by an expansion in terms of Chebyshev polynomials. Such a representation was introduced in Ref. [83] and has been used in many calculations since. Chebyshev expansions work very well with smooth functions. It is therefore advantageous to perform these expansions on a logarithmic grid for the logarithm of the function to be expanded. For a function f (x) this amounts to The t i are the Chebyshev polynomials and the R i /I i the respective coefficients.
As it turns out, our solutions are indeed smooth enough in the infrared and the ultraviolet momentum regions. However, in a short interval at intermediate momenta, we encounter large variations on rays in the second quadrant of the complex p 2 -plane. We deal with this situation numerically by splitting the radial distance from the origin on each ray into three intervals, [ 2 , x 1 ], [x 1 , x 2 ] and [x 2 , Λ 2 1 ]. Here, 2 is an infrared momentum cutoff and Λ 2 1 the ultraviolet momentum cutoff on each ray, already discussed above. The middle interval [x 1 , x 2 ] is bracketed close to the interval where the large variations occur. A set of boundaries that worked for the cases considered in this paper is given in table I. Furthermore, we optimized the number of Chebyshev polynomials N 1 , N 2 and N 3 needed for each interval; the resulting numbers can be found in the table as well. In total, we need to solve for the coefficients of 290 Chebyshev polynomials on each ray, which is a tremendous numerical task. For momenta with absolute value smaller than 2 , we use a constant extrapolation for the ghost dressing function G(p 2 ) and a constant extrapolation of the gluon propagator function D(p 2 ) = Z(p 2 )/p 2 into the deep infrared. Both extrapolations are well justified from the known behavior of the (decoupling case) dressing functions on the real axis and we assume this to work also in the complex plane. In fact, this can be checked by closely monitoring the behavior of the dressing functions for momenta larger than but close to 2 which was found to agree with the extrapolation.
FIG . 11. Two examples for the range of values for z (green line) for given values of x and y. The blue line represents the integration contour. The orange region is covered by the rays up to the one on which x lies. In the yellow region, the UV extrapolation is applied. When y is on the ray, all possible values for z lie also on the ray (and its UV extrapolation). When y is on the arc as in the plots, the required values for z form a line. The hatched region is only accessible once solutions for rays beyond the one on which x lies have been obtained. The lower half-plane can be accessed by complex conjugation of the dressing functions. The angle approximation corresponds to putting the green points all equal to y.
In the ultraviolet momentum region we used two types of extrapolation: first, one could use the analytic form of the known ultraviolet one-loop running of the dressing functions to extrapolate (see e.g. [53]), or, second, one could simply set the functions to a constant value from Λ 1 on. Both procedures lead to similar results and we settled with the simpler second option. We also need to integrate on the arc. This integration path generates momenta z in all regions of the complex momentum plane. To avoid this, we use the angular approximation Z(z), G(z) → Z(y), G(y) on the arc, which is known to work very well for ultraviolet momenta. For the main calculation, only the dressing functions on the rays are required. However, for some test calculations we also needed the dressing function between the rays. Since our rays are close to each other, we checked that linear interpolation between points on rays sharing a common distance to the origin works very well.

Numerical integration and iteration
For the numerical integration, some points require special attention which we discuss below. For the radial integration we use an ordinary Gauss-Legendre method separating the rays into four regions. For every external momentum x we use the three intervals of the Chebyshev expansion [ 2 , x 1 ], [x 1 , x 2 ] and [x 2 , Λ 2 1 ], detailed above. Additionally, we split the one region which contains x into two intervals. Thus, for small x, for instance, we have , with appropriate modifications for intermediate and large x. This is important to have many integration points close to the boundaries of the region and close to x ≈ y, where we pass through the only opening of the branch cut and encounter (depending on the angle) very small z in one of the denominators of the internal propagators. We typically choose 60 integration points in each interval. Another integration region is the interval [Λ 2 1 , Λ 2 ] on the arc. Since we use the angle approximation, this integration is not problematic. It turns out that the angular integral has to be treated with extra care, since it is related to (the generation of) cuts and branch points as explained in Sec. II. For any given pair of external and loop momenta x and y, respectively, the angular integral generates a potentially broad interval of values of z. To treat the evaluation of Z(z) and G(z) in a similar manner as with argument y, we perform the following procedure: For any pair of x and y, we monitor the interval of z tested by the angular integral. Whenever this interval crosses either x 1 or x 2 (the boundaries of the regions discussed above), we split the angular integral at these points. The total of 75 points used for the Gauss-Legendre integration of the angle are then distributed into these regions. Finally, we need to discuss details of the iteration procedure. The solution on the first ray/arc (only real and positive momenta) is obtained with standard techniques. We then use this solution as starting guess for the second ray/arc combination and iterate until convergence on the second ray is achieved. The solution for the second ray is then used as starting guess for the iteration on the third ray and so on. As explained in the main text, we renormalize the gluon DSE on the first ray/arc (only real and positive momenta); the corresponding value of Z 3 remains constant for all other rays. This procedure is not possible for the ghost propagator DSE, since G(s) at the infrared subtraction point s has a special meaning in connection with the family of decoupling solutions. To maintain connection to one particular member of the class of decoupling solutions dialed by G(s) on the first ray/arc, we use the Cauchy-Riemann condition discussed in section III B to determine the (complex) G(s) on each subsequent ray. Since, in particular in the infrared, the rays are very close to each other, the numerical error of this procedure is extremely small (we tested this explicitly on trial functions with various analytic structures in the complex plane).

Appendix B: Cauchy's integral formula and reconstruction
Cauchy's integral formula can be used to calculate the value of a holomorphic function inside a closed region via knowledge of the function on the boundaries: From this, one can directly derive the spectral representation of a propagator. Here we repeat this derivation but use different integration paths that correspond to the rays on which we calculate the propagators. This can be used as a test whether the numeric solution still respects analyticity.
In general, the analytic structure of a propagator can contain poles and a cut on the timelike axis. An integration contour of a circle at infinity thus needs to be deformed to take them into account: The integral at infinity (C ∞ ) vanishes and only the integral along the cut (C c ) and around the poles (C p ) remain. The latter leads to contributions of the residues of the n poles at z j : The first minus sign comes from integrating clockwise around the poles. The integral along the cut leads to This is the spectral representation for a propagator already shown in Eq. (9). The spectral density is defined as One can also change the contour such that it runs from infinity to the origin not along the timelike axis but along a ray at angle θ (and out at angle −θ), see Fig. 12. We consider the two contributions separately: If θ = π, we recover Eq. (9). Since we have a finite cutoff Λ 2 , we also add the contribution from the circle segment: The reconstructed propagator is then Appendix C: Three-gluon vertex models The different vertex models we employed in our studies are detailed here. For convenience, we define an auxiliary function δ = −9/44 is the anomalous dimension of the ghost propagator. The exponents of the dressings are taken such that this function behaves like (δ − γ), which is the anomalous dimension of the three-gluon vertex. The parameter a can be used to modify the IR behavior and was varied for testing. The baseline vertex is then defined as This is a Bose symmetrized version of the one introduced in Ref. [53] with an additional IR scale Λ 2 IR that serves to avoid the divergence for a = 3δ. The exponent accounts for the renormalization group improvement. We use as arguments x = p 2 , y = q 2 and z = (p + q) 2 as they actually appear in the gluon loop for external/internal momentum p/q. The factor 1/2 guarantees the correct behavior for large loop momentum. The integration kernel of the gluon loop contains the gluon dressing functions as Z(y)Z(z). To test the potential influence of and avoid possible problems from these terms, we discarded them in another vertex model and fixed the UV behavior by essentially modifying the exponent of the gluon dressing function from the original model: The potential poles induced by these arguments are irrelevant for the integration as they are on the negative real axis.
A final model we tested is the one from Ref. [53] which was also used for the previous calculation of the gluon and ghost propagators in the complex plane [26]. It is given by