Chaos in QCD? Gap equations and their fractal properties.

In the early 1980’s Benoit Mandelbrot pioneered the methodical study and computational visualization of the iteration of quadratic functions and began to cartograph the emerging fractal landscape [1] which subsequently has been named in his honor as the Mandelbrot set. With the advance of personal computers during the mid 80’s fractals gained broad attention, scientifically as well as in popular science.

In 1986, Leo Kadanoff, in an article with the title ”Fractals: Where’s the physics?”[2], expressed concerned curiosity about an understanding of fractal properties in physics which goes beyond the identification of fractal dimensions for certain problems. Kadanoff stated that without a better understanding of how physical mechanisms result in geometrical form it is difficult to trace types of questions with interesting answers. We wish to add that even in lack of such a deep understanding it is of course possible to find this kind of questions, as mentioned by Mandelbrot: ”I was asking questions which nobody else had asked before, because nobody else had actually looked at certain structures.” [3].

An example for this explorative approach is Hofstadter’s butterfly which is publicly less known. In 1976, ten years before Kadanoff asked his curious question and four years before Mandelbrot’s famous work on the quadratic map, Douglas Hofstadter observed what he called recursive structure in the computed spectrum of electrons in electromagnetic fields [4] which has been named after the visual appearance as Hofstadter’s butterfly. A first experimental confirmation of this theoretical prediction has been reported nearly twenty years later in 1997 [5].

There is no strict definition of what a fractal is but still most people would know one when they see it. Common descriptors of fractals refer to their non-analycity, self-similarity, non-linearity, iterative origin, chaotic behavior, non-integer (Hausdorff and other) dimension, to name a few. This paper has been motivated by the fact that QCD’s gap equations are by definition highly non-linear and self-consistent. Self-consistency equates the quantity of interest, our gaps, to a functional which depends on these gaps themselves. QCD’s gap equations are organized in a hierarchy of inter-dependencies of an infinite number of n-point Green-functions and it is at the heart of contemporary approaches in this field to identify methods which reduce this infinite number in a manageable way while preserving key features of QCD like dynamical mass generation and confinement. While one can argue how to obtain physically meaningful gap equations, viz. which set of approximations, truncations, etc. is the most reasonable, the self-consistent nature of these equations is not debated. Already at the seemingly simple level of 2-point Green functions for a single quark flavor appropriate truncation schemes allow to compute the mass spectrum of confined and deconfined quarks. The same methods allow for computation of meson and baryon spectra. Nothing of this is new and, although neither trivial nor brought to a final solution, it is in a structural sense reasonably well understood and dealt with in Dyson and Schwinger’s functional approach which proved to be a powerful tool to investigate the theories of QCD and QED. We refer to recent reviews for examples and more detailed information [6, 7, 8, 9, 10, 11, 12, 13].

Practitioners in the field of Dyson-Schwinger equations frequently deal with problems that can arise from their self-consistent nature. As an example, one technique to solve gap equations is by means of iteration starting from an initial guess. There is no guarantee for the convergence of such an iteration in general nor that the obtained solution is physical. In order to cover ’all possible’ solutions in this approach one would scan over different initial guesses. Typically, one can ’tame’ diverging iterations by damping the impact of the iteration itself. Instead of $g=F[g]$ one can write $g=\alpha g+(1-\alpha)F[g]$ where $g$ is the gap, $F$ a functional of the gap, and $\alpha$ a damping parameter close to but less than one, thus avoiding strong responses of $g$ to the iteration. One can wonder - we claim one should - whether it is justified to apply such an algorithm. It looks innocent in that sense that technically any solution of the original gap equation is a solution of the damped iteration equation. Nevertheless, at identical initial value both may give different answers and thus one can claim that the damping parameter might bear unwanted physical significance as it has been introduced ad hoc. What happens if the gap equation is allowed to iterate itself freely? We found only one, recently published, paper which asks exactly this question and comes to a clear conclusion: if the system is strongly coupled, chaos emerges and one can observe an infinite spectrum of ’unexpected’ gap solutions with increasing coupling strength[14]. In the paper we present, we give a brief explanation why these unexpected solutions actually should be expected. Further, we employ a model with momentum dependent gap solutions. In an iterative and inherently fractal context this led us on a surprising journey which answered not all but plenty of the questions we asked and at the end of which we are left to wonder whether looking at QCD as a fractal theory might be a key to understand confinement as an emergent fractal phenomenon.

Section II briefly motivates how iterative mapping generates new solutions of an equation while preserving the solutions of the non-iterated ’seed’ equation, Section III reviews the quark matter model by Munczek and Nemirovski (MN) in an extension for dense quark matter. We chose it for our exploration as it exhibits confinement and dynamical chiral symmetry breaking while being sufficiently simple to make it well suited for iterative mapping and analytic treatment. The following Section IV illustrates and cartographs chaotic features which emerge upon iteration of the gap equation, Section V is a cautious attempt to interpret physical meaning into the interplay of chaotic and non-chaotic structures we observe. Our study focuses at the structure of the mass pole. Appearance and disappearance of the mass pole are highly driven by chaotic behavior. Further, the mass gap itself is allowed to switch between different, usually distinct solutions. To our surprise, the physical properties of the iterative solutions provide a reasonable picture of how de-confinement could present itself in a model which possesses a gap equation with a single solution only. Finally, we estimate how a finite width gluon interaction could affect the observed behavior of the quark dispersion relation under iteration in Sec.VI before we conclude in Sec.VII.

We investigate possible consequences of chaos that appears in iterative solutions of non-linear and self-consistent equations in the complex domain. For clarity of what we consider physics and what is math we start with the latter and briefly review Mandelbrot’s fractal which is obtained by the iteration $z\xleftarrow{z_{0}=0;n\to\infty}f(z)$ with the explicit choice $f(z)=z^{2}+c$ to obtain the Mandelbrot set. We chose to use the symbol $\xleftarrow{z_{0};n}$ to have a distinguished notation for the iterative mapping process - specifying the number of iterations $n$ and the initial value $z_{0}$ - over the equal sign $=$ which appears in the analytic equation $z=z^{2}+c$. It is worthwile to look at the differences between these two. First, the polynomial equation has exactly two solutions $z_{1,2}$ for any given $c$ which are defined by the roots of the polynomial $P(z)=f(z)-z=z^{2}+c-z$. It is further easily seen that one can determine $c$ for a desired root $z_{0}$. For example, $P(z_{0}=0)=0$ if $c=0$.

In the iterative approach, each iteration generates a new polynomial,

etc. ad infinitum. There is one trivial but fascinating property of this infinite set of equations which essentially inspired the presented work. The left hand side of each of the previous equations has been set to $f_{i}(z,c)=z$ in order to obtain the next iteration $f_{i+1}(z,c)=z$. It is thus safe to state that the roots of $P_{1}(z,c)=f_{1}(z,c)-z$ are guaranteed to be roots of $P_{2}(z,c)=f_{2}(z,c)-z=f(f_{1}(z,c))+c$. As $P_{2}(z,c)$ is a 4th order polynomial, there are two more roots which, of course, did not appear for $P_{1}(z,c)$, a 2nd order polynomial. The important lesson to be learned is that for a self consistent non-linear equation $z=f(z)$ the iteration $z\xleftarrow{z;n}f(z)$ generates a new self consistent equation. While the solutions of the non-iterated equation remain a subset of solutions of the iterated equation, the iterated equations can develop additional solutions.

This is a peculiar, almost awkward situation, if one wishes to assign physical meaning to the original solutions of the equation $f(z)=z$. What makes these roots superior with respect to any of the iterative clones if all, original and clones, share these very same original solutions? Evidently, there is an infinite number of (iterated) functions which share the original roots. Is the original function with only these roots a superior or inferior function? Is it worth to ponder the meaning of the additional roots of iterated clones? Can we safely omit them? Do we miss important information when we ignore the duality of the gap equation as root defining equation and mapping rule? We decided to explore and ponder possible meaning.

One approach to move towards an understanding of QCD is based on evaluating QCD’s partition function by testing its response to external sources. This is the Dyson-Schwinger formalism which results in sets of coupled n-point Green functions. Out of these we are interested in the quark propagator which is obtained from the gap eguation

with the self-energy

Here, $m$ is the quark bare mass, $D_{\rho\sigma}(p-q;\mu)$ is the dressed gluon propagator and $\Gamma_{\sigma}^{a}(q,p;\mu)$ is the dressed quark-gluon vertex. This is the first of an infinite tower of gap equations which, without further approximations, couple back to this one. Further, there are similar equations for the dressed gluon-propagator and the quark-gluon vertex. Note that the gap equation is a self-consistent non-linear (in most cases integral) equation, $S^{-1}=F[S]$.

Within the Munczek-Nemirovsky model [15] the dressed quark-gluon vertex is approximated by the free quark-gluon vertex, $\Gamma_{\sigma}^{a}(q,p;\mu)=\frac{\lambda^{a}}{2}\gamma_{\sigma}$. Gap equations applying this approximation are referred to as rainbow gap equations. For the dressed gluon propagator the model is specified by the choice

Due to the $\delta$-function which in configuration space corresponds to a constant, this is a very simplified approximation of the gluon-propagator, specified by the coupling strength we set to $\eta=1.09$ GeV in accordance with [15]. For non-zero relative momentum $k$ the interaction strength in this model vanishes, thus making it super-asymptotically free. Furthermore, the infrared enhanced $\delta$-function is sufficient to provide for dynamical chiral symmetry breaking and confinement, both features of QCD which we wish to address. Finally, the $\delta$-function effectively turns the integral gap equation into an algebraic equation which can be solved analytically.

In order to obtain these solutions for the in-medium dressed-quark propagator one employs the general solution

Substitution into the dressed-quark gap-equation and appropriate tracing over the Dirac $\gamma$-matrices results in three coupled gap equations of which two (for $A$ and $C$) are identical:

We introduced $\tilde{p}^{2}=\vec{p}^{2}+(p_{4}+i\mu)^{2}$. In the chiral limit ($m=0$) one finds two distinct sets of solutions, one of them chirally symmetric and referred to as the Nambu phase,

whereas for the other solution, the Wigner phase, chiral symmetry is broken for $\mathcal{R}(\tilde{p}^{2})<\eta^{2}/4$,

For $\mathcal{R}(\tilde{p}^{2})>\eta^{2}/4$ the gap solution of the Wigner phase coincides with the Nambu solution. Note, that these solutions are obtained in Euclidean metric, but hold in Minkowski metric after a simple transformation, $\tilde{p}^{2}_{E}\to\tilde{p}^{2}_{M}$, with $\tilde{p}^{2}_{E}=\vec{p}^{2}+(p_{4}+i\mu)^{2}$ and $\tilde{p}^{2}_{M}=\vec{p}^{2}-(p_{4}+i\mu)^{2}$. Due to our interest in particle mass poles our investigation of the model is performed in Minkowski metric. For the next section however, the specific metric is not relevant; we will only work with the fact that $\tilde{p}^{2}$ is complex-valued and thus can be decomposed into real and imaginary part, viz. $\tilde{p}^{2}=z^{2}_{R}+iz^{2}_{I}$. We chose to label real and imaginary part with squared quantities as a reminder that they come in units of energy square.

Gap equations (5) and (6) lead to 4^th order polynomial equations with up to four distinct complex valued solutions at given $\tilde{p}^{2}$ for each gap. Generally, this is the whole solution space one would consider and the only task left is to identify the one physical solution. However, the self-consistent nature of (5) and (6) is evident and, according to our reasoning in the previous section there is a possibility for iterated functions with the same four and additional solutions.

Before we discuss our analysis, a few comments should be made. Defined by contact interaction in momentum space we chose a very simple model for the effective gluon propogator. For dressed-gluon propagators with finite width in momentum space the corresponding gap equations turn into integral equations. Thus, momenta couple and the simplicity of the MN model which we take advantage of for this exploration is lost. We address this issue in more detail in the last section of this paper.

As sketched in the Introduction, for models with sophisticated non-trivial interaction-kernel, iteration is a practical path to find gap solutions. We outlined before, that this leaves us with the possibility that iteration generates new functions which possess roots that correspond to solutions of the original gap equation and potentially an infinite number of additional roots.

We start our iteration from the non-interacting solution ($B_{0}=m$, $A_{0}=1$) and treat $\tilde{p}^{2}=z_{R}^{2}+iz_{I}^{2}$ as one would consider the constant $c$ for the Mandelbrot set $z\xleftarrow{}z^{2}+c$. For the moment, this reduces the number of independent variables from three ($\vec{p}^{2},p_{4},\mu$) to two. The result of such an iteration is shown in Fig.1 for the real part of the scalar gap $B$ at two different bare-quark masses of 10 MeV and 100 MeV, respectively.

Unlike the Mandelbrot fractal, this fractal does not diverge; chaos exhibits in domains in which the gaps for infinitesimally changes of energy and momentum take vastly different but finite values in a seemingly random pattern. This fractal region is contained within an almost perfectly shaped ellipsoid which we fit accordingly with

$(z_{R,0}^{2},R_{R}^{2},R_{I}^{2})$ differs slightly for $m=10$ MeV $(0.98,1.26,0.77)$ MeV² and for $m=100$ MeV $(0.99,1.28,0.80)$ MeV². The inner almond shape with less obvious chaotic behavior is well approximated by the same function with $(1.115,1.085,0.310)$ MeV² for $m=10$ MeV, and $(1.150,1.100,0.340)$ MeV² for $m=100$ MeV. As for the Mandelbrot set one would be ill advised to understand these figures as a valid representation of the fractal; the appearance of the fractal changes with each new iteration. We identify regions with identical periodicity, ranging from a stable, period one solution in the region outside of the covering ellipse over a period two region within the almond shape up to higher and higher periodicity in between these two regions. This is illustrated in Fig.3 for $m=100$ MeV for a periodicity of up to ten.

Keeping in mind, that the analytic gap equations possess four distinct solutions it seems interesting, that there is an extended stable domain (periodicity one) which favors one, and only one solution. We follow the gap solution along a vertical path at fixed $z_{R}^{2}$ and vary $z_{I}^{2}$. For reasons which become more clear at later stage we chose $z_{R}^{2}=0.1$ MeV². Along this line, one notices that one passes from an outer stable region into an inner stable region by traversing a small chaotic domain. This is illustrated in Fig.3. As within this chaotic domain the value of the gap function can change with each iteration, we plot all obtained values of $\mathcal{R}(B)$ over 300 iterations in gray scale according to how frequently a particular solution has been obtained. Evidently, there is a transition between two distinct analytic solutions of the non-iterated forth order polynomial gap equations. This result seems remarkable if one recalls how one would usually deal with different gap solutions for a given model: each solution is understood as a distinct phase, then one examines the stability of each individual solution and picks the energetically favored solution as the physical one. Upon iteration we are lead to a different conclusion. Although each of the analytic solution indeed is a solution of the gap equations, only one of them can be stable upon iteration at a given energy and momentum. However, the stable iterative solution over a finite range of energies can switch between distinct analytic solutions. It is further remarkable, that the iteratively stable solution is massive at low and massless at high energy. Amongst all the possibilities chaos seems to offer, this seems a very reasonable one. The notion of analytic solutions describing different phases however, is not supported from an iterative perspective; there is one, and only one, iterative solution to the gap equation.

Before we go into further interpretation of what this result implies, we wish to address a question related to the previous paragraph. Initially we remarked that our iteration starts from the non-interacting solution $A=1,B=m$. As we try to proceed as careful as possible, let us investigate the iterative stability of the four analytic gap solutions as plotted in the upper panels of Fig.4 where we show again the real part of the mass gap. The lower panel of Fig.4 shows the result after 300 iterations of these algebraic solutions as initial value. It is safe to say, that none of them is stable under iteration. Further, there is a visibly favored solution at large values of $z_{R}^{2}$ which does not depend on the initial gap that seeded the iteration. From a global perspective, the fractal keeps the general shape but shows differences for each different seed solution. This is to be expected and would happen in similar fashion to the Mandelbrot set if the initial value is arbitrarily changed.

Comparing the iterations to the algebraic solutions of the gap equations in the upper panel of Fig.4 one can graphically identify which of them is stable under iteration and in which domain. As seen, this is the case only for the positive mass-gap solutions 1 and 3 from Fig.4, as illustrated in Fig.5. In other words, although the chaotic domain will vary, the described features of Fig. 3 with respect to the analytic gap solutions do not critically depend on the chosen initial gap.

Up to this point we refrained from searching for meaning in our study. In spite of the fact that iterative mass gap solutions result in a large domain of chaotic behavior, which may or may not hide future surprises, we cannot help but wonder whether the switching between massive and mass-less gap solutions in the stable domains offers meaning. Before we go further, we want to recall that MN is considered to be a confining model. This is seen by the fact that the inverse propagator has no roots in the chirally broken phase and therefore integration over four-momentum does not pick up weight to generate a finite particle number. Hence, although confined quarks generate mass via chiral symmetry breaking, the absence of a mass pole results in the absence of a dispersion, viz. there is no explicit relation between specific momenta and energy. For the vacuum MN model, this is easily understood by the realization that in Minkowski metric $p^{2}+M^{2}(p^{2})$ has no real root if at any $p^{2}$, $M^{2}>-p^{2}=p_{4}^{2}-\vec{p}^{2}$. This running away of the mass in the chirally broken phase is exactly what happens in the MN model. However, as we have shown in the previous section iteration erases the distinction between chirally broken and restored phases and suggests that instead there might be a discontinuous gap solution which is confinement-like affected by dynamical chiral symmetry at small momenta, and at large momenta chirally unconfined-like and chirally restored. The transition between these domains is characterized by chaotic and unstable solutions (see Fig. 3).

At finite chemical potential the real poles of the propagator $A^{2}(\vec{p}^{2}-(p_{4}+i\mu)^{2})+B^{2}$ are represented by $\vec{p}^{2}-p_{4}^{2}+\mu^{2}+\Re(M^{2})=0$, with $M=B/A$. We note that the shift of the pole due to chemical potential should not be confused with the physical mass pole of the particle. This becomes evident if one considers an ideal non-interacting gas with $M$ being constant and real valued. For the purpose of this study we refer to the physical mass pole, defined by $\vec{p}^{2}-p_{4}^{2}+\Re(M^{2})=0$. From the definition $\tilde{p}^{2}=z_{R}^{2}+iz_{I}^{2}$ we identify the pole position in our contour plots as $z_{R}^{2}=\mu^{2}$ and $z_{I}^{2}=-2p_{4}\mu$ for an ideal particle with constant and real $M$. This represents a vertical line in our plots which does not depend on momentum and measures energy with increasing distance from the real axis. It shifts to higher $z_{R}^{2}$ with increasing chemical potential.

In Fig.6 we trace the physical mass pole in Minkowski metric by plotting the logarithm of the quantity $\left(p_{3}^{2}+M^{2}-p_{4}^{2}\right)^{2}$ which gives zero and hence a large negative logarithm at the physical mass pole. As the vertical axis does not depend on the mass ($z_{I}^{2}=2p_{4}\mu$) a vertical pole line indicates constant dressed quark masses. We observe the absence of such a well ordered pole structure within the chaotic domain. Since the vertical axis is a measure of the particle energy at fixed chemical potential one can conclude that the transition to the massive solution (Fig. 3) suppresses quasi-particle behavior in the infrared domain of the model. Again, we can trace the physical pole indicated by the vertical line and find $z_{R}^{2}=\mu^{2}-M^{2}$ since the pole is found at $p_{4}^{2}=p_{3}^{2}+M^{2}$. Following our elliptic fit of the outer boundary of the fractal domain this allows to determine the critical chemical potential where the infrared energy gap entirely disappears

We find $\mu_{C,IR}\approx 625MeV$ for $m=100MeV$ and $\mu_{C,IR}\approx 540MeV$ for $m=10MeV$. At these chemical potentials and beyond, quarks can be considered as completely chirally restored.

In order to estimate when mass-pole states can be occupied, we determine at which chemical potential the energy $p_{4}$ and the Fermi energy or chemical potential $\mu$ turn equal, that is when $z_{I}^{2}=\mu^{2}$ on the elliptic boundary of the fractal at the position of the physical mass-pole with $M=m$. We choose this scenario as this is the critical potential starting from where the particle energy is larger than the chemical potential and thus large enough to populate quasi-particle states. This is the case when

and holds for the light quark with $m=10MeV$ at $\mu_{m}\approx 359MeV$, for the heavy quark with $m=100MeV$ at $\mu_{m}\approx 432MeV$.

Although this is not a rigorous statement, one can roughly relate the critical chemical potential for the transition from a chirally broken into the restored phase to the in-vacuum dressed-quark mass. In our case, the situation is a bit different. We estimate a hypothetical chirally broken quark vacuum mass based on the previous estimate of the critical potential for the complete disappearance of the infrared gap by setting them approximately equal. Relating $\mu_{m}$ as onset of a deconfined, chirally restored quark phase with an estimate of the constituent quark mass seems to give rather reasonable results in comparison to other model calculations. This is interesting, considering that in the MN model the vacuum mass at zero 4-momentum is defined by the coupling strength $\eta$, which is of the order of 1 GeV.

It is noteworthy that our simple approach reproduces quantities related to the effective constituent masses at reasonable values. We state explicitly, that in this model constituent masses are nowhere realized for a physical particle, viz. an entity with a mass pole of that magnitude. We can compare the light quark critical chemical potential $\mu_{m}\approx 359MeV$ with the deconfinement critical potential obtained within the MN model in Euclidean metric with a value of $300$ MeV [16] or with subsequent work based on a widened version of the effective gluon propagator [17] which predicts deconfinement at a chemical potential of $380MeV$. There is a satisfying agreement of these values with ours. We point out though, that both of these models are defined within a different metric, slightly different bare quark masses and, most importantly, are based on entirely different assumptions. While the two previous papers employed distinct gap solutions and compare the pressure of the corresponding mass-less Wigner and massive Nambu phase, our approach results in only one gap solution which exhibits a transition from the Nambu to the Wigner phase through a chaotic domain as depicted in Fig. 3. Our quarks are either bare-mass quarks with poles or entities with a chaotic mass function or a dressed quark mass different from the bare mass with no associated pole. In the latter case there is a chaotic transition from dressed quark masses to bare quark masses with increasing energy.

We begin the final section of this paper with a plot of the particle pole in energy momentum space which we obtain by transforming $(z_{R}^{2},z_{I}^{2})$ to ($p_{3}=|\vec{p}|,p_{4}$) coordinates under explicit choices of the chemical potential as noted in Fig.7. Although this switch in representation does not provide additional information we find it instructive to provide an actual dispersion relation obtained from the iterative approach. In this example at a chemical potential of 700 MeV no chaotic behavior is visible and the dispersion is exactly that of a free quark at bare-mass 100 MeV. With decreasing chemical potential chaos emerges at energies higher than that of the expected (now absent) free particle dispersion. The actual dispersion branch is cut clean at some critical value (as we discussed in the previous section), thus illustrating our interpretation of the fractal boundary as the cause for a dynamical infrared cutoff below which quarks are mass-pole free.

We address a last question which relates to the fact that the MN model bases on the very particular choice of the effective gluon propagator as a $\delta$-function in 4-momentum space. The reason that we could easily perform the presented study bases on this Ansatz and the subsequent decoupling of momenta which allows to iterate point-wise for any given 4-momentum without coupling to other momenta. This might raise the suspicion, that momentum coupling could destroy the fractal structure we observed. In order to keep the simplicity of the gap equations but still get an idea about the stability of the emergent fractal we averaged at each point in our plane after each iteration step and thus mimicked some kind of momentum coupling. The averaging is based on Gaussian weights around a given point according to

where $w$ is the width of the Gaussian. To further simplify, we assumed that the widening only happens in the direction of momentum and energy, i.e. there is no widening perpendicular to the momentum.

As seen in Fig.7 the separation into chaotic pole-free, and non-chaotic mass-pole domains remains, even when we change the momentum dependence of the gluon from a delta function to a gaussian with a half width as much as 0.02 times the gluon mass. We find numerical evidence that this feature remains even with a width as much as 20% of the gluon mass, $\approx$ 200 MeV. This corresponds to a spatial width of about 1 fm; about the size of a proton. Based on this - certainly simplified - treatment of momentum coupling we conclude, that the statements we make in this paper may indeed survive a more complete treatment, involving self-interactions with globally coupled momenta - which has been our main concern prompting this final analysis.

As we have shown, a strictly iterative solution of the MN gap equations results in fractal gap structures which can be characterized by the existence of three qualitatively very different, yet co-existing domains of a single and unique gap solution: a bare-mass quark quasi-particle domain with physical mass-poles extending infinitely into the ultraviolet, a dressed-mass quark domain without mass-poles and hence no quasi-particle interpretation in the infrared, and a chaotic domain of transition between the first two phases. Remarkably, the two non-chaotic domains correspond to distinct analytic solutions which would usually represent individual phases with either dynamically broken or restored chiral symmetry. The fractal approach offers an alternative to this separation which is rooted in the iterative nature of the gap equation.

It is further noteworthy, that the iterative mass gap solution is always positive in the smooth, viz. non-chaotic domain of the fractal. The appearance of a chaotic boundary between two qualitatively different domains results in interesting properties:

I) The iterative approach provides an ultraviolet cut-off for the massive and mass pole free Nambu solution, as this solution appears only within the elliptic region of the ($z_{R}^{2},z_{I}^{2})$) plane. Thus the approach avoids the appearance of an infinitely increasing dressed-quark mass with increasing momentum and energy. In the MN model, this running mass results in the absence of mass poles for the massive gap solution and thus relates to confinement.

II) It provides an infrared cut-off for the bare-quark mass Nambu solution and thus ensures that quasi-particle states are not populated at small chemical potential although the quark can virtually exist as a quasi-particle with well defined dispersion.

III) Both cutoffs more or less coincide (as seen in Fig.5) although there is a transition region which is chaotic in nature. The resulting effective Nambu-UV/Wigner-IR cutoff depends dynamically on energy, momentum, bare mass, and chemical potential. As a side note we add that plotting gap solutions in the ($z_{R}^{2},z_{I}^{2}$) plane removes much of the dynamical arbitrariness and leaves the ratio of bare mass $m$ and coupling constant $\eta$ as the only ’true’ degree of freedom, viz. a change of the chemical potential $\mu$ would rescale the plot but cause no qualitative change, whereas plots like Fig.1 show indeed ’the’ gap solution at arbitrary chemical potential.

IV) At sufficiently large chemical potential bare-quark mass-pole states will form at energies which can be populated and thus physical quarks can exist as quasi-particle excitations.

A mechanism with these properties can be interpreted as a deconfinement mechanism. The appearance of one, and only one, iterative solution of the MN gap equations bears a certain elegance. First, by the very fact that there is only one gap solution with expected properties as the existence of only a positive mass gap, asymptotically restored chiral symmetry and the absence or appearance of physical mass poles. Next, it builds on distinct solutions which one would obtain in the non-iterative approach but gives them a new meaning by slicing them into a single new solution with the aforementioned properties.

A simple treatment of a widened, $\delta$-like gluon-interaction indicates that momentum coupling blurs chaotic domains but does not necessarily change the qualitative results we describe if the widening is moderate. As this study is an exploration and qualitative in nature we look forward to further analyses of this perspective on understanding confinement and the deconfinement transition as highly non-linear and to a certain extent possibly chaotic phenomena.

We are grateful to Prashanth Jaikumar, Pok Man Lo, and Craig D. Roberts for helpful comments and discussions which helped us a great deal to find order in chaos.