Super-Strong Coupling NJL model in arbitrary space-time dimensions

We study chiral symmetry breaking in the Nambu--Jona-Lasinio model regularized in proper-time in arbitrary space-time dimensions through an iterative procedure by writing the gap equation in the form of a discrete dynamical system with the coupling constant as the control parameter. Expectedly, we obtain the critical coupling for chiral symmetry breaking when a nontrivial solution bifurcates away from the trivial one and becomes an attractor. By increasing further the value of the coupling constant, we observe a second bifurcation where the dynamical solution is no longer an attractor, and observation that holds true in all space-time dimensions. In the super-strong coupling regime, the system becomes chaotic.

We study chiral symmetry breaking in the Nambu-Jona-Lasinio model regularized in propertime in arbitrary space-time dimensions through an iterative procedure by writing the gap equation in the form of a discrete dynamical system with the coupling constant as the control parameter. Expectedly, we obtain the critical coupling for chiral symmetry breaking when a nontrivial solution bifurcates away from the trivial one and becomes an attractor. By increasing further the value coupling constant, we observe a second bifurcation where the dynamical solution is no longer an attractor, and observation that holds true in all space-time dimensions. In the super-strong coupling regime, the system becomes chaotic. Quantum chromodynamics (QCD) [1,2] is nowadays accepted to be the theory that describes strong interactions among quarks and gluons. It possesses two important and opposite features, asymptotic freedom [1] in the high energy domain, where quarks move freely at short distances. Perturbative QCD is the appropriate framework to approach this regime. On the other hand, there exists quark confinement [3] in the low-energy and largedistances regime. In this case, the fundamental degrees of freedom bind in hadronic bound-states (mesons, baryons and exotic), the coupling becomes strong and the theory highly non-linear. A number of non-perturbative frameworks have been developed to study low-energy QCD, among which we find Lattice QCD [4], Schwinger-Dyson equations (SDEs) [5,6] (see Refs. [7][8][9][10][11][12][13] for reviews of SDEs in hadron physics) and other field-theoretical approaches, like the functional renormalization group [14][15][16], as well as effective models. Nambu-Jona-Lasinio (NJL) model [17,18] is a favorite one to study the properties of hadrons mostly in connection with chiral symmetry breaking (see, for instance, the reviews in Refs. [19,20]), and it is the general topic of this article. The model is non-renormalizable and thus, there are several schemes in which the regulator and coupling are selected to match static properties of pions (see, for example, Ref. [21]). Here, we select to regulate the gap equation in the proper-time (PT) scheme. Momentum integrals become Gaussian and thus are readily performed; the regulator comes in the PT integration, as opposed to the more frequently employed hard momentum cutoff regularization scheme [22] (the UV regulator in momentum space in the latter case is identified with the IR cut-off in proper-time integrals). Furthermore, by introducing an IR cut-off in momentum (UV in PT), we mimic confinement in the model [23][24][25]; the second cut-off removes quark-antiquark production thresholds, in the deep infra-red region of the quark propagator. This procedure has an added advantage that the resulting quark propagator naturally fulfills the Axial-Vector Ward identities for bound states. Furthermore, the PT regularization scheme allows an easy incorpora-tion of plasma effects, including finite temperature, density and external magnetic fields. Among other things, these extreme conditions make the effective dynamics of quarks reduced from the usual 4-dimensions in vacuum; a strong magnetic field reduces two dimensions, whereas a heat bath at high temperature reduces one dimension. On top of that, the coupling of the original theory effectively changes when quarks interact in a medium in such a manner that these external agents, temperature and magnetic field, compete with each other to either inhibit or promote the dynamical breaking of chiral symmetry, namely, the effective coupling is dressed by the medium in such a way that its strength can be enhanced or diluted. Therefore, it is natural to ask for the behavior of the solution to the gap equation in several space-time dimensions and for a broad range of values of the coupling (see, for instance, Refs. [19,[26][27][28][29][30][31][32][33]) and explore how scenario of ChSB changes. This is the goal of the present article.
Regarding its parametric structure, NJL model has a single coupling constant that must exceed a critical value G c to describe chiral symmetry breaking. That is to say, for G > G c , a non-trivial solution to the quark gap equation M = 0 bifurcates from the trivial one M = 0. In the chiral limit (see below), the gap equation can be cast in the form where the form of f (M ; G, Λ, . . .) depends on the regularization scheme through the cut-off(s) Λ, whereas the coupling constant G plays the role of a control parameter in the language of discrete dynamical systems [34] and the dots stand for any other possible parameters under consideration, like temperature, density, magnetic fields and so on. A solution M = 0 turns out to be an attractor in the sense that regardless the initial condition M 0 , after a finite number of iterations, the succession M 0 , M 1 = f (M 0 ), . . . , M n = f (n) (M 0 ) is such that M n lays within a arbitrarily small vicinity around the solution M . Thus, this procedure is useful to find the solution to the gap equation with a desired precision merely by considering more and more iterations. In Ref. [35], for example, by using a hard momentum cut-off regularization scheme, the authors cast the gap equation for the NJL model in terms of a discrete dynamical system with the coupling of the model as control parameter and interpret each iteration as a self-energy insertion in the quark propagator. They identify the critical value of the coupling for chiral symmetry breaking, G c , consistent with other approaches and find the dynamically generated mass to be an attractor even for large values of the coupling constant G > G c . These findings were generalized by some of us in Ref. [36] for different regularization schemes. In that work, we observed, not surprisingly, that for the three-dimensional (3D) and fourdimensional (4D) hard momentum cut-off regularization schemes, the procedure converges to the dynamical mass in a finite number of iterations and this solution continues to be an attractor even in the super-strong coupling regime. On the other hand, for the Pauli-Villars (PV) and PT regularization schemes, we noticed that a second bifurcation appears for a coupling larger than G c such that the solution ceases to be and attractor and a double-period orbit develops. For even larger values of the coupling, the system becomes chaotic. This means that due to the non-renormalizability of the model, care has to be taken when resumming self-energy corrections and is the point we pursue to stress in the manuscript: When high-momentum modes in the self-energy are not removed to start with, the resummation is no longer stable. Our goal in this article thus becomes to explore in a closer look the circumstances under which the first and second bifurcations of the solutions to the gap equation in the super-strong coupling regime take place in arbitrary space-time dimensions. For that purpose, we have organized the rest of the manuscript as follows: In Sect. I we present the basics of the NJL model in vacuum for arbitrary space-time dimensions. The general iterative procedure to solve the gap equation is discussed in Sect. II. Numerical findings of the evolution of the first and second bifurcations in arbitrary space-time dimensions are discussed too. The analysis for the confining variant of the model is discussed in Sect. III. We present a summary of our findings and perspectives in Sect. IV.

I. NAMBU-JONA-LASINIO MODEL
The Nambu-Jona-Lasinio (NJL) model [17,18] was first introduced to explain chiral symmetry breaking for nucleons inspired by the phenomenon of spontaneous symmetry breaking in superconductors even before the QCD era. In modern considerations, it is now regarded as an effective model of interactions among quarks and gluons in which the one gluon exchange between the quarks is replaced by four-fermion contact interaction, as shown in pictorial form in Fig. 1. The Lagrangian of the NJL model for two light quark flavors N f = 2, and three colors N c = 3, is given by [19,20] where the first term is the matter (Dirac) part of the Lagrangian with m = diag(m u , m d ) denoting the current light quark mass, which in the chiral limit we set m = 0. The second term is the four-fermion interaction with (ψψ) 2 is the scalar piece of the interaction, while (ψiγ 5 τ ψ) 2 is axial-vector part. Here τ are the Pauli matrices acting in isospin space and G denotes the coupling constant. The above Lagrangian (1) respects all the global symmetries of the the full QCD Lagrangian in chiral limit [19,20]. Chiral symmetry can be explicitly broken via a finite current quark mass. It can also be broken spontaneously. In this regard, a starting point to study the latter scenario is the gap equation, shown pictorially in Fig. 2. From a Schwinger-Dyson equations (SDEs) point of view, we start from the most general form of the dressed quark propagator where Σ(p) represents the quark self-energy, which in QCD corresponds to the expression with D µν (q 2 ) denoting the gluon propagator, g 2 the QCD coupling constant, Γ ν (k, p) is the dressed quark-gluon vertex, λ a are the SU (3) Gell-Mann matrices. Mapping one-gluon exchange diagrams to a contact interaction in the NJL model equivalent to consider in Eq. (3). Then, it is straightforward to verify that there are no wavefunction renormalization effects in this model, and that the quark mass function is merely a constant such that the dressed quark propagator takes the form with M denoting the momentum independent dynamical quark mass, which is self-consistently determined from the gap equation The symbol Λ stresses that the model is not renormalizable and we need to regulate the integrals involved. Within the PT regularization scheme, we exponentiate the denominator of the integrand in (6) such that where in the last step we have performed the Gaussian momentum integrations. Notice that Λ and λ are regulators with mass dimension -1. Thus, the gap equation in arbitrary dimensions can be written as and the effective coupling With these ingredients, and for the sake of exploring the parameter dependence of the gap equation, we rewrite Eq. (8) in dimensionless form as with µ = M Λ, µ 0 = mΛ and f (µ) = E D 2 (µ 2 ). Before proceeding to analyze Eq. (11), let us consider a variant of the regularization in Eq. (7) that has been recently introduced and has the form [23][24][25] This double cut-off PT regularization has the advantage of mimicking confinement. Unlike the standard regularization in Eq. (7), the form of the quark propagator regularized as in Eq. (12) is consistent with the Axial-Vector-Ward identities when used in the Bethe-Salpeter equation [24,25]. This model is dubbed Contact Interaction (CI) model, and its coupling in four-dimensions is parameterized as G CI = 4πα ir /m 2 G , where m G = 800 MeV is a dynamical gluon mass parameter [38] and α ir = 0.93π is the interaction strength. For this CI model, the gap equation reads which in terms of the dimensionless variables µ and µ 0 described above and defining ρ = λ/Λ, we re-write as with G CI D = DG CI /(2) D−1 (π) D/2 . In the following, we study the solutions of Eqs. (11) and (14) through an iterative procedure.

II. SOLUTION OF THE GAP EQUATION
A. Analysis of NJL model in D = 4 In the following, and for the sake of illustration, we consider the solutions to the gap Eq. (11) in the chiral limit, µ 0 = 0, in the case D = 4. The ideas presented here shall be extended to other cases below. Let us commence, inspired by the analysis of the logistic map [34], by observing that the strength of the coupling, G 4 in our notation, dictates the number of solutions and its nature in the dynamical systems sense. In Fig. 3 we plot the r.h.s. of the gap equation as a function of µ for three values of the coupling in Eq. (10). The upper panel corresponds to G 4 = 0.8, G 4 = 1 and G 4 = 1.3. We also plot the curve y = µ as the l.h.s. of the gap equation and look for the intersection of both these curves at fixed G 4 . We observe that for G 4 < 1, the only intersection on the corresponding two curves is at the origin µ = 0 and the r.h.s curve always lies below the diagonal line. By observing the cobweb plot (also shown), we notice that the trivial solution corresponds to a fixed point. It is an attractor in the sense that given an arbitrary initial guess µ i , in a finite number of iterations, we can find the iterations in a vicinity of the origin. The r.h.s. curve is taller as G 4 grows bigger, being tangent to the diagonal line when G 4 = 1 near the origin. As the coupling increases, the r.h.s. curve intersects in two points, the origin and a new point µ f = 0 which turns out to be fixed and an attractor, as seen in the cobweb construction. By still demanding the coupling to be larger (lower panel in 3 we show the cases G 4 = 3, G 4 = 3.5 and G 4 = 4), the point µ f ceases to be fixed at the point µ * when the slope of the r.h.s. curve reaches its critical value f (µ s ) = −1, which turns out to be the case for G 4 = 3.5, resulting in a flip bifurcation [34]. This is a similar test performed in the logistic map to describe precisely this type of bifurcation. For larger couplings (G 4 = 4 in Fig. 3) the cobweb plots shows a 2-cycle orbit, observed as the corners of the black square. A doubling-period cascade develops for larger values of G 4 until the system becomes chaotic, as observed in Fig.4. This picture is observed in arbitrary dimensions (see Fig. 5) and also with the CI model (see Sec. III) for different values of the second cut-off. Nevertheless, we would like to emphasize that the bifurcation map of Fig. 4 does not correspond to the actual solution of the gap equation. As mentioned before, the procedure of iterating the gap equation is equivalent to the loop expansion of the propagator and with each iteration, the number of loops added to the propagator increases by 1. Following this idea, at one hand, the convergence of the iterations for relatively small couplings means that higher loops do not contribute and thus the iterative procedure converges to the actual solution to the gap equation in this case, the plot corresponds to the solution of the gap equation. In Fig. 4, it is shown as a stable branch before the flip bifurcation develops. On the other hand, when the coupling increases, the non convergence of the iterations points out that every added loop contributes and thus we got the bifurcation diagram even though for the gap equation itself there is only one solution, which, by the way, never exceeds the natural scale of the model, namely, the cut-off. Below we present a detailed analysis in these cases.

B. Analysis in arbitrary dimensions
The bifurcation map described above replicates in arbitrary dimensions, as can be observed in Fig. 5. Therefore, an analysis similar to the logistic map can be done to identify the first and second bifurcations in the map. To this end, we start from the gap eq. (11) and take the derivative with respect to µ in both sides, namely The critical coupling for the first bifurcation hence appears when we evaluate the above expression at µ = 0, when the nontrivial solution bifurcates from the trivial one. Thus, the critical coupling for chiral symmetry breaking in arbitrary dimensions is given as and it is plotted in Fig. 6. We observe that the slope of the curve of this critical coupling as D increases is uniform, suggesting a linear dependence with D.
A standard criterion to locate flip bifurcations in dynamical systems such as the logistic map is found when the slope of the map under study reaches the critical value -1 [34]. We adapt this criterion which in our case consists in locating the position µ * where the function The variation of µ * in different dimensions is shown in Fig. 7. It monotonically increases with D till it saturates at large number of dimensions. The value of the critical coupling for the second bifurcation, G * D , is We plot this coupling in Fig. 6. It linearly rises with the dimension D, but changes slope around D 4. For dimensions D lower that 4, the stable branch, namely the portion of the plot where the iterative procedure converges to the actual solution of the gap equation before the second (flip) bifurcation appears, is shorter than for D > 4, as can be seen from the separation between the two couplings. This means that the larger the dimension, the fixed point branch associated with chiral symmetry breaking is more stable.

III. CONFINING MODEL
We can extend the above analysis to the confining version of the model described in Eq. (14) in the chiral limit µ 0 = 0. In this case, the derivative of the gap equation should be taken for f CI (µ). First, we notice as the cutoff ρ grows, the critical value of the coupling for the first branching, namely saturates at a fixed D for its value as in the ordinary model, eq. (11). This is shown in Fig. 8, where we observe that G 0 CI D → G 0 D as ρ → ∞. At fixed ρ, the variation of G 0 CI D with the dimension is shown in Fig. 9. A linear growth is evident to the eye.
The critical coupling for the second bifurcation, G * CI . We have fixed ρ = 1.1. is, as before, found by locating the value of µ * at which in the form and is also shown in Fig. 8. Again, it grows linearly in two regimes separated around D = 4, where the slope changes such that the first branch of fixed-point solutions is stable for a larger interval of values of G CI D . The scale µ * is shown as a function of D in Fig. 10 for fixed ρ. We notice that the impact of the cut-off ρ is less important for large D, and that the variation of µ * is less for a smaller cut-off.
Finally, we notice the appearance of doubling period orbits for larger values of the coupling. In Fig. 11 we show this is the case for fixed ρ in different arbitrary dimensions. We observe that the effect of the dimension influences directly the length of the dynamical chiral symmetry branch, between the first and second bifurcations. It also can be seen in the maximum and minimum values of the orbits at fixed G D in the chaotic domainbifurcations are wider for low dimensions in the range of couplings plotted. In Fig. 12 we depict a similar bifurcation pattern for fixed D = 4 and different values of ρ. In this case we notice that for larger ρ, orbits saturate and tend to the same pattern as in the ordinary NJL model shown in Fig. 4. All these observations are consistent with the statement that the chaotic behavior is a robust feature of the PT regularization of the model, with one or two cut-offs in arbitrary space-time dimensions.

IV. SUMMARY AND PERSPECTIVES
By using an iterative procedure for the solution of the gap equation for the NJL model, we explore the behavior of the dynamical quark mass as a function of the coupling strength in two variants of the PT regularization of the model. We naturally find a critical value G 0 D for the coupling required for breaking chiral symmetry in arbitrary dimensions. Such critical value specifies the minimal value of the coupling where the dynamical mass bifurcates away from the trivial value µ = 0. The dynamical solution is a fixed point and behaves as an attractor for the iterative procedure as expected for an interval of values of G. Nevertheless, a second bifurcation occurs when the coupling exceed the value G * D , where the dynamical solution ceases to be an attractor for the iterations. Larger values of the coupling trigger a cascade of doubling-period orbits [34] until the gap equation becomes chaotic in the dynamical systems sense. This happens for arbitrary dimensions and even avoiding quark production thresholds in the regularization.
On physical grounds, external agents like temperature or density dilute the strength of the interaction and induce a dimensional reduction of the system when the temperature and/or chemical potential are very large. A uniform magnetic field, on the other hand, plays an opposite role and in spite of the fact of inducing a dimensional reduction, it increases the effective strength of the coupling and hence it could drive the effective strength of G to very large values. Therefore, the simple minded exercise carried out in this article allows to address the qualitative behavior of the dynamical mass as a function of the effective coupling in different arbitrary dimensions. The stable branch of the bifurcation has a finite length in this view and thus, the appearance of the second bifurcation sets bounds on the validity of the PT regularization in the sense of stability of the iterative procedure to solve the gap equation.
At this moment, it remains unclear whether the interpretation of every iteration as a self energy insertion in the propagator advocated in [35] is physically reliable or if it is regularization scheme dependent, if the gap equation within the PT (and PV and, perhaps, others) regularization schemes do not admit an iterative solution for super-strong couplings or if the nature of the chaotic behavior found in this article has implications in actual physical systems. Research along these lines is being carried out and results will be presented elsewhere.

V. ACKNOWLEDGMENTS
AA acknowledges the Department of Physics Gomal University (Pakistan) and IFM-UMSNH, (Mexico) for support. AA also thanks to the staff of IFM-UMSNH, (Mexico) and Paola Rioseco for the hospitality during the visit where part of this work was carried out. AM and AR acknowledge Consejo Nacional de Ciencia y Tecnología (Mexico) for support under Grant 256494.