Bifurcations in the RG-Flow of QCD

Bifurcation analysis is used to study an effective model of QCD4 with four-fermi interactions. Our analysis supports the scenario of a fixed point merger at the lower edge of the conformal window. This indicates square root scaling of the anomalous scaling dimensions of the fermion fields just above the lower edge and exponential scaling just below. We also predict existence of new fixed points in this model whose (dis)appearance may indicate transitions of the flow within the conformal window. Furthermore, we make new predictions for the critical value (Nf/Nc)crit at the lower edge. We also obtain exotic spiraling flows that are generated by complex scaling dimensions of the effective four-fermi interactions. Finally, we extend the model by adding a scalar field that couples with a Yukawa interaction term and study the modifications it causes to the RG-flow.


Introduction
Renormalization group describes the flow of a theory in the space spanned by the coupling constants called the "theory space". In quantum field theory renormalization group flows are governed by a set of beta functions that result in ordinary differential equations. The flow induced by these beta functions provide information about the change of coupling constants and anomalous field dimensions with the energy scale. An interesting aspect of the RG flows are the fixed points where the beta functions vanish. In particular one is interested in determining the positions of these fixed points and their stability properties.
In addition one would like to obtain possible renormalized trajectories generated by the JHEP07(2019)075 beta functions. Since the RG flows are generated by ordinary differential equations, it is natural to study them with the tools developed for the theory of dynamical systems. In this article, we use analytical and numerical 1 methods for bifurcation analysis to study bifurcations in renormalization group flows, as recently proposed by Gukov [3]. Bifurcation theory, which we summarize in appendix A, categorizes topological changes in a dynamical system and is therefore an ideal tool to study the appearance and disappearance of fixed points. Using this categorization, we can not only determine easily when fixed points appear and disappear on the RG flow, but also obtain information about the nature of these transitions where the fixed points (dis)appear. As discussed in [4], there are three typical ways in which fixed points of an RG-flow (dis)appear: • Collision of a fixed point with non-trivial coupling constant and a fixed point with a trivial coupling constant. At this collision the non-trivial coupling changes sign and the stability of both fixed points changes.
• Annihilation of two fixed points with non-trivial coupling constants. The annihilating fixed points have different stability types.
• A divergent fixed-point as one or multiple coupling constants run off to infinity.
The first and the second scenarios are examples of standard bifurcations. The first is called a transcritical bifurcation, while the second is referred to as the saddle-node bifurcation. A mathematical definition of and the conditions for which these bifurcation occur, can be found in appendix A. For the first scenario it is essential that one of the fixed points has at least one trivial coupling constant. It then collides with another fixed point for which the same coupling constant is non-trivial. The non-trivial fixed point changes sign in this coupling constant at the bifurcation. In some cases a negative coupling constant means that the fixed point moves into an unphysical region, as is for example the case for the gauge coupling in QCD. For this reason this scenario was described in [4] as a merger of a trivial and non-trivial fixed point, where the non-trivial fixed point disappears. Furthermore, this scenario is likely to occur in systems of beta functions, where one of the beta functions has a trivial root. An example is the beta function of the gauge coupling in QCD, that always has a root for trivial gauge coupling.
The first two scenarios also determine the scaling of the anomalous dimensions of operators that cross marginality at the bifurcation. It is easy to show that the first scenario always induces a linear scaling of the anomalous dimension at the fixed point [3], while the second scenario implies square root [3] and exponential (also known as BKT or Miransky scaling) scaling of the anomalous dimensions [4] at the two sides of the merger.
Finally, the third scenario above with diverging fixed points is much harder to describe using techniques from bifurcation theory, since no local bifurcation occurs. However, this scenario is mostly encountered in a perturbative analysis of the RG flows which becomes invalid at the divergence. It turns out that this last scenario will in some cases be replaced JHEP07(2019)075 by the second one, if the perturbative methods are improved for example through the use of an effective field theory approach. Furthermore, this scenario may be ill-defined in nonperturbative approaches, as the CFT data becomes ill-defined. It's thus unlikely to happen in any physical system. 2 In this paper, we apply bifurcation analysis to an effective model of QCD as suggested in [3]. We consider QCD 4 with gauge group SU(N c ) with an arbitrary number of colors N c and flavors, N f , as the bifurcation parameters. The Lagrangian is given by where G A are the gluon field strengths and ψ i are the fermions of flavor i. Furthermore a, b, c are used to label the color in the fundamental representation, and A, B, C is used for the color in the adjoint representation. We have We will consider the massless case, m i = 0 ∀i, where the Lagrangian becomes invariant under a global SU(N f ) L × SU(N f ) R × U(1) symmetry which can be broken to SU(N f ) V by formation of a chiral condensate. The two-loop beta function for the gauge coupling is given by [5] β [2] g ≡ Λ with α g := g 2 (4π) 2 and are the scheme independent beta-function coefficients. This beta function has a (degenerate) trivial root at α g = 0, and a non-trivial root at α g = − b 0 b 1 . The non-trivial fixed point becomes negative and therefore unphysical for 5) and where the latter bound runs between 3.4 and 2.6 for 1 < N c < ∞, thus it is always below the first bound for any N c . The region between these bounds, where a non-trivial fixed point 3 exists, is known as the conformal window (cf. figure 1), and its existence was first JHEP07(2019)075 discussed in [5,6]. Within the conformal window the non-trivial fixed point is attractive in the IR limit, while the trivial fixed point is attractive in the UV. Above the conformal window, the trivial fixed point becomes an IR attractor. Since α nt g becomes small close to x = 11/2 higher order contributions to the beta function (1.2) can be neglected, hence the value of upper bound is robust. However, close to the lower edge α g becomes very large, invalidates the above estimate and (1.6) becomes unreliable. In particular, other fixed points resulting from higher order terms might influence this critical value.
The type of phase transition that occurs at the lower edge of the conformal window, that we call x crit , and its value has been a central problem in the literature. Use of Schwinger-Dyson equations [17,19,21,22] generically suggest that the theory is in the confining phase with chiral symmetry breaking for x below x crit , see figure 1. As for the value of x crit , estimates for N c = 3 are typically in the range 6 x crit 12 [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23], and follow from the various techniques (cf. table 6 of [3] for a recent overview of various studies and their prediction). A more recent study [24,25] uses Potts models to discuss the behaviour at the lower edge of the conformal window. The transition at the lower edge could be first order [26,27] or infinite order [4,22], and can be analyzed through several methods. This is an important open problem, as characterization of this phase transition will lead to a better understanding of the strong interactions and might be relevant for model building beyond the standard model [27]. Furthermore, functional RG and holographic approaches indicate that the conformal window is continuously connected to the quark gluon plasma phase at finite temperature [28][29][30]. Therefore studies of the conformal window may provide crucial information on the quark gluon plasma and vice versa.
Several works in the literature [9,24,25,27,28] point toward the possibility that the IR stable fixed point α g = −b 0 /b 1 annihilates with another fixed point of the RG flow at x crit in the form of a saddle-node bifurcation (second option in the list above). This would in turn indicate that a UV irrelevant operator, such as a 4-fermi interaction, crosses marginality precisely at x crit [3] as already anticipated in the early literature [20,21,23]. In this paper, we investigate this possibility by considering an effective model for QCD 4 that includes four-fermi interactions and determine all fixed points within this model in and around the conformal window. Furthermore, we discuss how they appear and disappear as the bifurcation parameters N c and N f are varied. Finally we study the effects of adding a scalar operator with a Yukawa type interaction to our results in this effective theory.

JHEP07(2019)075
The paper is organized as follows. In the next section we introduce the effective model with the four-fermi interactions, derive the beta function equations both for finite values of N c and N f and in the Veneziano limit, and study the fixed point structure using bifurcation theory in detail. In section 3, we extend the model including the scalar operator. We conclude and discuss our results and list possible future work in the last section. Appendix A provides a short introduction to bifurcation theory in dynamical systems and appendix B and C detail the calculation of the beta-functions.
2 An effective model for QCD

Beta functions in an effective model for QCD
We study the theory in the Wilsonian picture. For this, we write down the effective action, which is given by the Legendre transform of the Lagrangian (1.1) with m i = 0. We impose the symmetries obeyed by the Lagrangian on the effective action and expand the effective action in the quark fields ψ ai up to the 4-quark interactions, as was done in [15]. This yields where L 4f denote the four-fermi interactions and we redefined the gauge fields One can write down four independent four-fermi operators [31][32][33][34] respecting the symmetries. We choose a basis as in [15,34]: where η = γ ψ is the anomalous dimension of the fermion field and One can calculate the beta functions of the four-fermi interactions using the Wetterich equation given in appendix B.3. The beta function for the gauge coupling α g cannot be derived from the effective action (2.1) easily as one needs an expansion in the gluon fields in the effective action [15] which complicates the calculation substantially. Instead, we can start from the perturbative equation (1.2) and add the effective four-fermi interactions to JHEP07(2019)075 it. We provide this calculation in appendix B.2 and obtain the same result as in [15]. The complete set of beta-functions is then given by where we defined the rescaled variables (2.5) The same model was studied in [15]. However, the calculations performed in appendix B leading to (2.4) yield slightly different beta functions from the ones presented in [15]. The beta functions presented in this paper agree with those in [15] in the Veneziano limit, that was only analyzed in [15]. In this paper we'll also consider the finite N regime, and we thus need the full set of beta functions.
Let us further define and perform the rescaling

JHEP07(2019)075
which makes the set (2.4) amenable to the Veneziano limit where N → ∞ and x remains finite. After the rescaling we obtain Before analyzing this model, let us shortly discuss its validity: one can obtain an exact model by either using perturbative methods and calculate corrections to all orders or one can use an effective field theory approach, which contains infinitely many effective operators. In this paper, we apply a mix of the two methods as was done in [15]. The effective interactions will then incorporate parts of the higher loop correction to the gauge beta function. As can be found in appendix B.2, we have constructed the beta functions in such a way that there is no overlap between the perturbative and effective corrections to the beta functions.
The beta functions for the effective interactions are calculated by the Wetterich equation and thus exact at one loop. However, errors are introduced due to the truncation of the effective action. For simplicity we have only considered the effective four Fermi interactions in this model. Inclusion of running couplings for higher dimensional effective operators could destabilize fixed points found in this model or introduce new fixed points and bifurcations.
However, as long as we stay within the class of models of the same dimension (number of beta functions) that have the same symmetry, we can guarantee that the results will persist for all sufficiently small model deformations. More precisely, if the corresponding bifurcation is non-degenerate, then any perturbed beta function will have it for nearby parameter values. This is also true if we add higher-order terms, but then new equilibria and their bifurcations can appear for large amplitudes.
Since model deformations generally scale with the values of the coupling constants, our results are perturbatively stable as long as coupling constants are small. We will take a bound on our parameter values of α i < 1, |g i | < 1, and notice in our results when this bound is exceeded. Our results should not be trusted beyond this bound. Below this bound the one-loop effects will dominate over the higher-loop effects, but higher-loop effects might yield significant contributions to the beta functions, which will be ignored in JHEP07(2019)075 this toy model. Finally, for α i 1, |g i | 1, higher loop contributions will be negligible. Furthermore, the results of our model can be generalized to models including higher order effective operators, as long as the coupling constants for these operators are small and taken at their fixed point value.

Veneziano limit
In this section, we study the Veneziano limit of the model (2.7) as was done in [3] and [15].
where we have defined t = − ln(Λ/Λ 0 ). It is important to note that all higher order contributions in the perturbative theory are either suppressed or give constant contributions in the Veneziano limit with the rescaled variables. We now observe that the first three equations decouple from the last two. Therefore, we can reduce analysis to the system (2.9) The domains of the variables and the parameters are given by x ∈ R + , α ∈ R + 0 and g S , g V ∈ R. The beta function for α g has a double root for α t g = 0 and another root for . Therefore, the fixed points of the RG-flow lie on these manifolds. The manifold defined by α g = 0 is an invariant set of the system, sinceα g vanishes for α g = 0.
The behavior of the fixed points projected on the (x, α g ), (x, g S ), (x, g V ) planes and (x, α g , g S ), (x, α g , g V ) and (x, g S , g V ) spaces are shown in figure 2. We use the following color coding in this figure: • Solid red line: stable node, 3 negative eigenvalues (IR attractor).    . The fixed points of the RG-flow of (2.9) projected on the (x, α g )-plane, (x, g S )plane, (x, g V )-plane, (x, α g , g S )-space, (x, α g , g V )-space and the (x, g S , g V )-space. LP: Limit Point (Saddle-Node Bifurcation). BP: Branching Point (Transcritical Bifurcation). The equilibria on the α g = 0 manifold all have a trivial eigenvalue, and therefore do not fit the color coding as described above. In order to use the color coding, we define the sign of a trivial eigenvalue on this manifold by approaching the equilibrium along a line of constant g S and g V from positive but small α g . The trivial eigenvalue will then approach 0 from either positive or negative values. If the trivial eigenvalue is slightly positive for 0 < α g 1, we define the trivial eigenvalue to be positive and vice versa. The transcritical bifurcation at x = 1 is located in the subspace α g = 0, while the other transcritical bifurcations are located on the intersection of the subspaces α g = 0 and α g = α nt g . In the latter case, both a trivial and a non-trivial fixed point 6 are involved. At all the non-trivial fixed points, the fixed point has a positive eigenvalue, when α g < 0, that becomes negative when α g > 0, while the trivial fixed point has a negative eigenvalue that becomes positive. At the branching point at x = 5.5, a stable node with α g > 0 is created. This indicates the upper edge of the conformal window, or creation of the Banks-Zaks fixed point, as discussed in the introduction of this section. Finally, at x = 4.05 two non-trivial fixed points collide and disappear through a saddle-node bifurcation (limit point). This saddle-node bifurcation could indicate the lower edge of the conformal window through a fixed point merger, as was discussed in [15] and [3]. The two other non-trivial fixed points diverge for x 2.5, close to the divergence of the large N c limit in equation (1.6).

JHEP07(2019)075
A few characteristic orbits of the RG-flow in the Veneziano limit are visualized in figures 3a, 3b and 3c.
In figure 3a, a phase portrait on the invariant plane α g = 0 is shown at x = 5. We find 4 fixed points, and 4 red critical heteroclinic orbits, connecting those fixed points, and marking the boundaries of an invariant set. All orbits inside this set have a finite UV and finite IR fixed point, while orbits outside the invariant set have a diverging IR limit or diverging UV limit. The structure of this flow does not change, when x is varied in the region x > 1. Only the location of the fixed points and connecting orbits will change JHEP07(2019)075 (a) Phase portrait of the RG-flow in the (g S , g V )-plane for α g = 0 in the Veneziano limit at x = 5.
(b) Phase portrait of the RG-flow in the (α g , g S , g V )-space in the Veneziano limit at x = 5.
(c) Phase portrait of the RG-flow in the (α g , g S , g V )-space in the Veneziano limit at x = 4. continuously. At x = 1, the two fixed points with smallest value of g V collide and exchange stability, via a transcritical bifurcation, yielding a similar phase portrait.
In figure 3b, a phase portrait of the flow in the 3-dimensional space is shown at x = 5. Here, all 8 fixed points are shown, together with a few heteroclinic orbits connecting them. The red heteroclinic orbits form a skeleton of a 3-dimensional invariant set. This invariant set is bound by six invariant planes. The solid red lines lie on the intersections of these planes. The dashed red lines indicate the presence of other intersections, of which the exact locations have not been found numerically. All orbits inside this set have a finite UV and IR fixed point, while orbits outside the invariant set have a diverging IR limit, a diverging UV limit or both. More precisely, the bulk of the invariant set consists of structurally stable orbits going from the fixed point at (α g , g S , g V ) = (0, 0, −0.333) to the fixed point at (α g , g S , g V ) = (0.032, 0.003, 0.000), and the boundary planes consist of non-stable heteroclinic orbits between the other fixed points. Within the invariant set the theory is in a chirally symmetric phase, while outside the set the chiral symmetry is broken. The structure of this flow does not change when x is varied in the region 4.05 < x < 5.5. Only the location of the fixed points and connecting orbits will change continuously. At x = 4.05 the two fixed points shown at (α g , g S , g V ) = (0.032, 0.003, 0.000) and (α g , g S , g V ) = (0.073, 0.322, −0.084) disappear through a saddle-node bifurcation.
In figure 3c, a phase portrait in the whole space is shown at x = 4. Here, only 6 fixed points are left. The six fixed points together span two unstable bounded 2-dimensional invariant sets. The solid red lines border these sets, and the dashed red line indicates the location, where the last border is expected, which has not been found numerically. The 3-dimensional invariant set, which was present at x = 5, has been broken down at the saddle-node bifurcation to two 2-dimensional invariant planes. The structure of the flow does not change when x is varied in the region 0 < x < 4.05, but the two non-trivial fixed points diverge when x is lowered towards x = 2.5. The orbits outside of the invariant planes all diverge in the IR limit. The eigenvector with 0 eigenvalue at the saddle-node bifurcation 7 is mainly directed in the g S direction, as can be found in table 2. Therefore, the divergence is the strongest in the g S direction, indicating that the chiral symmetry is broken due to the scalar part of the four-fermi interaction, as discussed in [35]. Furthermore, we notice that there still exist chirally symmetric trajectories with a zero UV and finite IR limit, but these are unstable, and therefore need a very precise fine-tuning within this model. In the model out of the Veneziano limit there are 2 bifurcation parameters N c and x = N f /N c , which allows for more interesting types of bifurcations. Again, we are interested in the renormalized trajectories of the theory and the invariant sets enclosed by the renormalized trajectories, indicating the existence of chirally symmetric theories. We have numerically found that the model contains up to 20 different real fixed points with α g ≥ 0 for integer values of (N c , N f ) ∈ [1, 10] × [1, 100]. The beta function for α g yields two roots α g = 0, α g = α nt g := N 2 where the trivial root has multiplicity two. The (α g = 0)-hyperplane is an invariant set of the model, which allows to study the model without gauge interactions separately. Setting
(b) Saddle-node bifurcations in the (α g = α nt g ≥ 0)-subspace. Black lines: saddlenode bifurcation, that was also found in the Veneziano limit. α g = 0, we find up to 12 fixed points, which can be continued using MatCont. Two of those disappear through a saddle-node bifurcation at x > 7.5, which is above the conformal window. The 10 remaining fixed points play an important role for the RG-flow for α g = 0. For every of these ten fixed points, we find a related fixed point on the manifold α g = α nt g . For large values of x all the non-trivial fixed points have α g < 0, but for smaller values of x they have α g > 0, and become physically relevant. On the intersection of the two hyperplanes, where α g = α nt g = 0, i.e. when x = 11 2+3g V , all non-trivial fixed points collide with a trivial fixed point, through a transcritical bifurcation. The lines at which these transcritical bifurcations occur are shown in figure 4a. Using the transcritical bifurcations, we relate all the non-trivial fixed points to a trivial fixed point.
The 10 non-trivial fixed points with positive α g exist for large values of x. However, for small values, they either disappear through a saddle-node bifurcation or diverge in one or multiple couplings. In figure 4b the saddle-node bifurcation lines are shown. Here, we see that two pairs (black lines) of fixed points will disappear through a saddle-node bifurcation which is connected to the saddle-node bifurcation that was found in the Veneziano limit. 8 The two bifurcation curves are close together for large N c , but split up at small N c . One of the curves starts to increase steeply around N c = 2, while the other starts decreasing around N c = 6, and diverges in the coupling constants at N c = 2. Therefore the saddlenode bifurcation of this pair does not exist for N c = 1, and the fixed points diverge instead. Then there is a third pair of fixed points (upper blue line) that disappears through a saddlenode bifurcation for all N c around x = 4.8. However this value starts to increase around N c = 4. A fourth pair disappears through a saddle-node bifurcation when 1 < N c < 16, JHEP07(2019)075 diverges for N c ≥ 16, while at N c = 1 the fixed points will not exist for positive α g . The last two fixed points do not disappear through a saddle-node bifurcation, but will diverge in multiple couplings at small x.
Using these saddle-node bifurcation curves, we can partition the 10 non-trivial fixed points into 5 pairs. Since all 10 non-trivial fixed points are related to a trivial fixed point through a transcritical bifurcation curve, we can partition the 20 fixed points into groups of 4 consisting of 2 trivial and 2 non-trivial fixed points. It turns out that the heteroclinic orbits between the fixed points within each of these sets have a similar topology, which is as shown in figure 5: at first, for large x, the two non-trivial fixed points lie below the α g = 0 hyperplane, and the α g = 0 hyperplane is locally attractive in the IR limit, then, at a specific value of x, one of the non-trivial and one of the trivial fixed points collide through a transcritical bifurcation and the non-trivial fixed point gets a positive α g creating a locally attractive IR limit set bounded by orbits between the 3 fixed points with α g ≥ 0. Next, at smaller specific x, the other non-trivial fixed point will collide with the other trivial fixed point and also get positive α g value. After this, at small x, the non-trivial fixed points will either disappear through a saddle-node bifurcation or diverge, making the α g = 0 hyperplane locally repulsive in the IR limit.
We can also visualize the structure of the flow on the invariant (α g = 0)-space. A graphical representation of the fixed points and a few critical orbits is shown in figure 6 for (N c , N f ) = (3,15) and (N c , N f ) = (3, 3). We notice that two points (labeled by 5 and 7) do exist for (N c , N f ) = (3, 15), but do not exist for (N c , N f ) = (3,3). This is due to a saddle-node bifurcation in the (α g = 0)-space at (N c , N f ) ≈ (3, 5.9). Furthermore two fixed points (labeled by 6 and 9) have exchanged stability through a transcritical bifurcation at (N c , N f ) ≈ (3, 7.7). The curves corresponding to these two bifurcations are shown in figure 7e and 7a respectively. In the remainder of this section, we'll use the numbers as in figure 6 to label the fixed points. Since we have related every trivial fixed point uniquely to a non-trivial fixed point, we give the non-trivial fixed point the same number as its trivial counterpart. However, we add the label 'a' to the trivial fixed points and the label 'b' to the non-trivial fixed points.
We have numerically found all local bifurcations 9 involving these 20 fixed points for α g ≥ 0, and these are shown in figure 7. This figure includes the saddle-node and trans-JHEP07(2019)075 In addition, the non-trivial fixed points at N c = 3 are shown as a function of x for α g > 0 in figure 8. In this figure the solid red line indicates the Banks-Zaks fixed point labeled with number 10b. The branching points are the transcritical bifurcations, where the non-trivial fixed points branch of their trivial counterpart. Al these branching points are points on the solid red transcritical bifurcation curves in figure 7. Furthermore, the limit points represent the saddle-node bifurcations at N c = 3 and lie on the solid blue lines in figure 7.
We use the bifurcation diagram in figure 7 to distinguish and visualize several phases within the RG-flow generated by our model (2.4). For this, we set N c = 3, and take N f ∈ {11, 12, 13, 15}. The topology on the invariant (α g = 0)-space will be the same for these four values of x, since there are no bifurcations in this region. Only the location of the fixed points will change continuously. Projections of the fixed points and a few heteroclinic orbits in this 4-dimensional invariant subspace, at α g = 0, are shown in figure 9. The fixed points and a few heteroclinic orbits in the α g = α nt g subspace are shown in figures 10, 12, 13 and 14. In figures 11 and 15 these figures are shown in a graphical representation.    ) and α g = α nt g . Solid line: heteroclinic connections. Dashed lines: heteroclinic connections that are expected to exist, but have not been found numerically. Left: projection on the (g S , g V , g V2 )-space. Right: projection on the (g V , g V1 , g V2 )-space.
At N f = 15, we find 10 fixed points and a few heteroclinic orbits connecting them, both in the α g = 0 (figure 9) and the α g = α nt g hyperplane (figure 10). The red lines shown in the figures form part of a skeleton of a closed set in the 4-dimensional subspaces α g = 0 and α g = α nt g . These two 4-dimensional sets on α g = 0 and α g = α nt g form a 4dimensional boundary of a larger 5-dimensional invariant subspace of the full 5-dimensional space. The other 4-dimensional subspaces that bound the full 5-dimensional invariant set are spanned by orbits going from α g = 0 towards α g = α nt g . Dashed red lines indicate parts of the skeleton that have not been found numerically, but are expected to exist. The flow within the 5-dimensional invariant set is everywhere directed from the α g = 0 towards the α g = α nt g > 0 subspace. This set consists of structurally stable renormalized trajectories that have finite UV (with trivial α g ) and IR limits (with strictly positive α g ) indicating the theory to be in a symmetric phase. Outside the set, orbits have diverging IR and/or UV limit, and theories will have a broken chiral symmetry. The boundaries consist of unstable orbits between other fixed points, and therefore indicate unstable symmetric theories.
(b) Topology of the RG-flow at α g = α nt g corresponding to figure 10. When the flavor number is reduced just below N f = 15, at N f = 14.7 and N f = 14.6, two pairs of fixed points disappear through a saddle-node bifurcation. Because of these bifurcations, the invariant set of orbits indicating symmetric theories is reduced. However, the most attractive IR fixed point with dim(W s ) = 5 still exists ( figure 12). This indicates that the invariant set is locally attractive in the IR limit. Reducing the flavor number even further to N f = 12 ( figure 13) and N f = 11 (figure 14) more fixed points disappear through a saddle-node bifurcations. The invariant set of orbits at α g = α nt g is reduced at every step, and therefore the set of asymptotically stable orbits with finite UV and IR limit becomes smaller, making the set of possible symmetric theories smaller. For N f ≤ 12.2, there is no fixed point left with dim(W s ) = 5. Therefore, the remaining invariant set is not an IR attractor anymore. Furthermore, the dimension of the invariant set with finite trajectories reduces, when more fixed points have disappeared through a saddle-node bifurcation. When the invariant set splits up into multiple invariant sets with dimension smaller than 5, symmetric theories only exist within this model under extreme fine-tuning. This is a strong indication that the lower edge of the conformal window has been crossed.
In figure 9, there is another peculiarity visible. Here we see a heteroclinic orbit with spiraling behavior indicating the existence of renormalized trajectories with complex scaling behavior in and close to the α g = 0 hyperplane. Equilibrium 5 and 7 are associated to this behavior.

Operators crossing marginality
We have found 4 saddle-node bifurcation curves (cf. figure 4b, 7) with non-trivial value of α g . All four indicate the disappearance of two fixed points when the value of x is JHEP07(2019)075       Table 3. Eigenvalues and critical eigenvectors crossing marginality of the saddle-node bifurcation with non-trivial α g at N c = 3.
decreased. Along these curves we can identify the eigenvector corresponding to a zero eigenvalue as the relevant operator that crosses marginality. These critical eigenvectors at N c = 3 are reported in table 3. We have found that these operators are very dependent on x, but not on N c . From the table we conclude that 3 of the saddle-node bifurcations have 2 complex eigenvalues at N c = 3, while one has only real eigenvalues. We notice that the normalized eigenvectors corresponding to the complex eigenvalues are all directed along the four-fermi couplings, and have very small contributions (|η αg | < 0.007) along thê e αg direction, meaning that this complex scaling behavior involves operators, which are linear combinations of the effective four-fermi interactions. Furthermore, we see that the operator crossing marginality is a linear combination of all operators. However the scalar four-fermi interaction is most relevant in three of the saddle-node bifurcations, while in the saddle-node bifurcation of equilibria 5b and 7b the V 1 interaction is most relevant.

Degenerate bifurcations
We have encountered nine bifurcation points with higher codimension. An overview can be found in table 4. Six of these points lie in the physical uninteresting region N c ∈ (0, 1), while three of them (CP, SNT1 and IP1) may have physical relevance, since N c ≥ 1. Those points will be discussed below. The point labeled with 'CP' corresponds to a cusp bifurcation and the points labeled with 'SNT' correspond to saddle-node-transcritical bifurcations, as discussed in appendix A. The point 'TCT' labels an intersection of two independent transcritical bifurcation curves and corresponds to the normal form ẋ = αx + x 2 , y = βy + y 2 . (2.10) The points 'IP' label more degenerate bifurcation points. At the cusp bifurcation 2 saddle-node bifurcation curves with non-trivial value of α g meet. The first curve is associated to the pair of equilibria (6b,8b) and the second to (6b,9b).
The point SNT1, describes a saddle-node-transcritical bifurcation, of a saddle-node curve that lies in the subspace α g = 0 and involves the equilibria (2a,3a), and two transcritical curves that also lie in the α g = 0 subspace and involve the equilibria (2a,4a) and (3a,4a). IP1 is a highly degenerate bifurcation point. The model reduced on the invariant subspace α g = 0 has a pitchfork bifurcation at this point. However, in the 5-dimensional model the branches of the equilibrium curves meeting at this pitchfork bifurcation are all transcritical bifurcation curves of the equilibria (6a,6b), (8a,8b), (9a,9b) and (6a,9b). In addition to this pitchfork bifurcation, we find a saddle-node-transcritical bifurcation at this point. The curves associated to this saddle-node-transcritical bifurcation all lie in the α g = 0 subspace. The saddle-node curve involves the equilibria (8a,9a) and the branches of the transcritical bifurcation involve (6a,8a) and (6a,9a). On top of this, there is a saddlenode bifurcation curve with non-trivial α g intersecting at this point. This curve involves the equilibria (6b,8b). Together with the branches of transcritical bifurcation curves of JHEP07(2019)075 (6a,6b) and (8a,8b), that were contained in the pitchfork bifurcation, we could also view this point as saddle-node-transcritical bifurcation of these points. In short this point is highly degenerate and can be unfolded into two saddle-node-transcritical bifurcation and a pitchfork bifurcation of transcritical bifurcation curves.

Discussion of the results
We distinguish 10 pairs of fixed points, that may be relevant for the physics in and below the conformal window. These points can be divided into three different sets: {1, 2, 3, 4}, {5, 7} and {6, 8, 9, 10}. Here, we used the numbering as in figure 6. The first set contains the most repulsive fixed points, while the last set contains the most attractive fixed points.
The set of pairs of fixed points {6, 8, 9, 10} consists of points that were also found in the Veneziano limit. The non-trivial points of the pairs 8 and 10 disappear through a saddlenode bifurcation at x ≈ 4 for large values of N c like in the Veneziano limit. Furthermore, the non-trivial points of the pairs 6 and 9 diverge for N c ≥ 16, and the trivial points of the pairs 6 and 9 have a transcritical bifurcation at x = 1 for large values of N c . This behavior was also found in the Veneziano limit.
Fixed point 10 is the most attractive and the point with non-trivial α g in this pair describes the IR limit of the structurally stable renormalized trajectories that exist within the invariant set. Furthermore, it is continuously connected to the Banks-Zaks fixed point that was found in the Veneziano limit. We see that like in the Veneziano limit the non-trivial IR point disappears through a saddle-node bifurcation with the non-trivial fixed point of . This saddle-node bifurcation is a good candidate for the lower edge of the conformal window in this model. This prediction is higher than, but comparable to the predictions summarized in [3]. Furthermore, it is close to a recent prediction on a similar model including effective interactions [34]. In addition, we see that the transcritical bifurcation of the fixed points of pair 10, where the most attractive fixed point obtains a positive α g is constant at the ratio x = 5.5. This line indicates the upper edge of the conformal window, and the same value was found in the Veneziano limit. This is consistent with earlier predictions, as was expected, since α g is small, so the two-loop perturbative beta function of α g without four-fermi interactions should give a good approximation. The resulting estimate for the conformal window is shown in figure 16a.
Furthermore, we see that the two fixed points that were diverging in the Veneziano limit for small values of x will disappear through a saddle-node bifurcation for values 1 < N c < 16. This saddle-node bifurcation takes place for ratios below the conformal window for N c > 5, and for 1 < N c ≤ 5 it takes place in or above the conformal window, while at N c = 1 these fixed points do not exist for strictly positive α g . This saddle-node bifurcation curve affects the closed invariant set of structurally stable orbits, and could therefore be relevant for the physical behavior of QCD in the conformal window.
For N c ≥ 16, we expect similar behavior as in the conformal window, where these two fixed points diverge around x ≈ 2.5. Here, the fixed points form extremal points of a 5-dimensional closed invariant set, which is discontinuously broken down to a smaller invariant set, when the saddle-node bifurcation of the pairs 8 and 10 takes place.

JHEP07(2019)075
On the other hand, for 5 < N c < 16, the points no longer diverge, but also disappear through a saddle-node bifurcation. However, this happens for a value of x which is below the lower edge of the conformal window. Therefore, the behavior within the conformal window will not be changed much by this saddle-node bifurcation. In the region 3 ≤ N c ≤ 5, this saddle-node bifurcation takes place within the conformal window, resulting in a discontinuous change in size of the closed invariant set of structurally stable orbits (or renormalized trajectories with finite IR and UV limit) at the bifurcation. I.e. it does not affect the range of the conformal window, but it induces an discontinuous change in the set of symmetric theories, which could indicate a transition within the conformal window in the sense that a non-empty set of finite renormalized trajectories, which are present near the upper edge of the conformal window, are not present close to the lower edge and have diverging IR limits instead.
At N c = 2, the saddle-node bifurcation takes place above the conformal window, and therefore this discontinuous transition will not take place within the conformal window. However, the set of symmetric theories in the conformal window will be smaller than at larger values of N c . At N c = 1, this saddle-node bifurcation only takes place for negative α g , but not for positive α g .
The set of pairs of fixed points {1, 2, 3, 4} contains the four pairs of fixed points that are the most repulsive of the 10, and are therefore most relevant for the UV physics. Pair 1 is the most repulsive pair and therefore describes the UV limit of the renormalized trajectories within the invariant set of structurally stable orbits. At a curve x > 8 the fixed points within this pair undergo a transcritical bifurcation, as can be seen in figure 7c. Below this curve the most repulsive of the pair is the trivial fixed point making the theory asymptotically free. Above this curve the α g = 0 hyperplane is repulsive, and therefore the theory always diverges in the UV limit.
The non-trivial fixed points of the pairs 1 and 2 both diverge when x is decreased. Furthermore, we see that two of the non-trivial fixed points, the ones from the pairs 3 and 4, disappear through a saddle-node bifurcation for N c ≥ 2 this bifurcation happens for smaller values of x than the saddle-node bifurcation which we took as the lower edge of the conformal window in figure 16a. We note that this bifurcation curve also converges to x ≈ 4 in the Veneziano limit, where the saddle-node bifurcation was found in the Veneziano limit.
We find that both this curve and the saddle-node bifurcation curve of the non-trivial fixed points 8 and 10 correspond to the saddle-node bifurcation that was found in the Veneziano limit, since the beta functions for g V 1 and g V 2 have two fixed points at this point, which shows that there are actually 2 saddle-node bifurcations in the Veneziano limit, that differ in the values g V 1 and g V 2 .
The saddle-node bifurcation curve of the non-trivial fixed points of the pairs 3 and 4 is another candidate for the lower edge of the conformal window (figure 16b). If we take the saddle-node bifurcation curve of points 3 and 4 as the lower edge of the conformal window, the saddle-node bifurcation curve would indicate a discontinuous change in the invariant set of structurally stable orbits indicating the symmetric theories. Furthermore, this would imply that the lower edge of the conformal window is just below N f = 12 at N c = 3, as was found in lattice studies such as [11], and from Schwinger-Dyson equations with ladder JHEP07(2019)075  resummations in [19][20][21][22][23], where the lower edge was predicted at Finally, we notice that there is a last set of two pairs of fixed points, {5, 7}, which was not found in the Veneziano limit. These fixed points have complex eigenvalues with relatively large complex part compared to the real part. This induces a spiraling behavior in the system. These points disappear through a saddle-node bifurcation both on the α g = α nt g and the α g = 0 subspace. These fixed points might have a physical interpretation or might be a consequence of the approximate nature of our model and would then be destabilized further when higher dimensional effective operators are taken into account.
Furthermore, the other fixed points (except for the pairs 8 and 10) can have complex eigenvalues, but for those the complex part is smaller than the real part, making the spiraling behavior less visible. Also, the eigenvectors corresponding to the complex eigenvalues are largely directed into the direction of the effective four-fermi couplings. The complex eigenvalues indicate that when the orbits move away (towards) the fixed point there is both a scaling, given by the real part of the eigenvalue and a rotation, given by the complex part of the eigenvalue. The scaling determines the velocity of flowing away from (towards) the fixed point, while the rotation indicates that the relevant operator, which is a linear combination of the basis operators rotates in theory space. In our case this rotation takes place in the subspace spanned by effective operators, since the operators with complex eigenvectors are mostly directed along the effective operators.

JHEP07(2019)075
In conclusion, we find two saddle-node bifurcations that could indicate the lower edge of the conformal window, in the sense that both involve the merger of a non-trivial fixed that is created through a transcritical bifurcation at x = 5.5, and thus corresponds to the pertubative Banks-Zaks fixed point. Both scenarios are shown in figure 16. These two fixed points merge with other non-trivial fixed points, which are created at values x > 5.5 through transcritical bifurcations. The two possible lower edges of the conformal window are then given by (N c , N f ) = (3.00, 12.2) and (N c , N f ) = (3.00, 11.3). Apart from these 2 pairs non-trivial fixed points, we find 3 more pairs of non-trivial fixed points that are all created through transcritical bifurcations at x > 5.5 and either disappear through a saddle node bifurcation or diverge. 10 Furthermore, we find that most bifurcations occur at almost constant x(N ), for values of N ≥ 3.

QCD 4 with a scalar field
In this section we extend our previous model with a meson-like scalar field Φ that couples to the fermion fields with the coupling constant y. The field Φ is invariant under SU(N c ) and transforms in the adjoint representation of SU This extension is based on an extension mentioned in [34]. The Lagrangian is extended with a massless scalar field and a Yukawa coupling with the fermion fields yielding extra terms in the Lagrangian: The set of beta functions is then extended to a set of 6 differential equations in 6 variables and 2 parameters, which is calculated in appendix C:

JHEP07(2019)075
where we have defined α y = y 2 (4π) 2 . We rescale the couplings as in (2.6), together with the rescaling N c α y → α y , and we define the ratio x := N f Nc , and N := N c . We find

The Veneziano limit
We take the Veneziano limit (N → ∞), which yields (3.4) The first 4 equations decouple from the last two and we can reduce analysis to (3.5)

JHEP07(2019)075
The domain of the variables and parameters is given by x ∈ R + , α g , α y ∈ R + 0 and g S , g V ∈ R. The beta function for α g has a double root for α t g = 0 and another root for . Therefore, the fixed points of the RG-flow lie on these manifolds. The manifold defined by α g = 0 is an invariant set of the system, sinceα g vanishes for α g = 0.
In addition, the beta function for α y has a root for α t y = 0 and another root for α nt y = 3αg+2g S 2+x . Therefore, the fixed points of the RG-flow lie on these manifolds. The manifold defined by α y = 0 is an invariant set of the system, sinceα y vanishes for α y = 0.
The location of the fixed points projected on the (x, α g )-plane, (x, α y )-plane, (x, g S )plane and the (x, g V )-plane is shown in figure 17, and the projections on the (x, α g , α y )space, (x, α g , g S )-space, (x, α g , g V )-space, (x, α y , g S )-space, (x, α g , g V )-space and the (x, g S , g V )-space in figure 18. Here, we use the following color coding: • Solid red line: stable node, 4 negative eigenvalues (IR attractor). The equilibria on the α g = 0 manifold all have a trivial eigenvalue, and therefore do not fit the color coding as described above. In order to use the color coding, we define the sign of a trivial eigenvalue on this manifold by approaching the equilibrium along a line of constant α y , g S and g V from positive but small α g . The trivial eigenvalue will then approach 0 from either positive or negative values. If the trivial eigenvalue is slightly positive for 0 < α g 1, we define the trivial eigenvalue to be positive and vice versa. 11 Since β αy = 0 if α y = 0, we find the same fixed points and bifurcations as in the previous model without scalar field (cf. figure 2). In addition, we now find five fixed points with non-trivial value of α y . We find two additional branching points and one additional limit point, as can be seen in table 5. Furthermore, at 3 transcritical bifurcations, that were found in the α y = 0 model, we now find three equilibria that intersect and exchange stability, corresponding to the normal form ẋ = αx + x 2 , y = αy + y 2 . (3.6) One of the additional fixed points branches from the transcritical bifurcation at x = 1.000 that was already found in the model without a scalar, and diverges rapidly in α y . In addition, one new fixed point branches from one of the equilibria with non-trivial value of α g at x = 3.926, which also diverges rapidly in α y . Furthermore, at the branching point at x = 8.173, an extra fixed point branches of with a non-trivial value of α g and α y , JHEP07(2019)075  . The equilibria of the RG-flow of (3.5) projected on the (x, α g , α y )-space, (x, α g , g S )space, (x, α g , g V )-space, (x, α y , g S )-space, (x, α y , g V )-space and the (x, g S , g V )-space. LP: Limit Point. BP: Branching Point.  Table 5. Bifurcations found in the Veneziano limit of QCD 4 with a scalar field (3.5).

JHEP07(2019)075
which diverges when x is lowered towards x ≈ 4. Another new fixed point with non-trivial α y exists in the scalar model. For this point α y → 0 if x → ∞ and α y → ∞ when x is lowered towards x ≈ 2. From this point another fixed point branches at x = 11.647. This point disappears trough a saddle-node bifurcation at x = 5.497. The other point that disappears at this saddle-node bifurcation branches at x = 5.5 from the trivial fixed point (α g , α y , g S , g V ) = (0, 0, 0, 0). This point is the stable node in the system. This stable node disappears through a saddle-node bifurcation at 5.497. We therefore find a conformal window with non-trivial for α y at x ∈ (5.497, 5.501), which is smaller than a prediction on a model without effective four-fermi interactions, where x ∈ (5.24, 5.5) [34].

Small N c regime
In this section, we analyze the complete model (3.3) for small N c . The beta function for α g has two roots: where the trivial root has multiplicity 2. The trivial root makes the α g = 0 hyperplane an invariant set of the model (3.3). The beta function for α y also has two roots: The trivial root makes the α y = 0 hyperplane an invariant set of the model (3.3). Since α y = 0 is an invariant hyperplane, we find the same bifurcations as in previous section. However, on top of those, there exist fixed points and bifurcations in the α y = 0 space. We find up to six fixed points with α g , α y = 0 and up to six fixed points with α g = 0, α y = 0, depending on the values of (N c , x). As in the previous section, we find that fixed points with non-trivial α g always collide on the intersection α nt g = 0 with a fixed point with trivial JHEP07(2019)075 Figure 19. Topology of the RG-flow between 4 fixed points.
value α g = 0 through a transcritical bifurcation. This happens at where the second condition is a reduction to the model without scalar field. Similarly, fixed points with non-trivial α y always collide on the intersection α nt y = 0 with a fixed point with trivial value α y = 0 through a transcritical bifurcation. This happens at where only the positive root is a physically relevant solution. In this way, one can relate most fixed points with (α g > 0, α y > 0) to a fixed point with (α g = 0, α y > 0), a fixed point with (α g > 0, α y = 0) and a fixed point with (α g = 0, α y = 0). We find that the orbits between these points are always as indicated in figure 19. Furthermore, a non-trivial fixed point 12 with positive α g and α y can appear through a transcritical bifurcation at (α g = α nt g = 0, α y > 0), at (α g > 0, α y = α nt y = 0) or at (α g = α nt g = 0, α y = α nt y = 0). These transcritical bifurcation curves are shown in figure 20 by dash-dotted red lines, dash-dotted black lines and solid red lines respectively. When x decreases non-trivial fixed points will either diverge in one or multiple couplings or disappear through a saddle-node bifurcation with another fixed point. We find three saddle-node bifurcation curves with (α g = 0, α y = 0) as shown in figure 20a, 20b and 20c. They're shown together with the related 13 saddle-node bifurcation curves that were found for α y = 0 in the previous section. In addition, a saddle-node bifurcation curve with (α g = 0, α y = 0) related to the saddle-node bifurcation in figure 20c is shown. 12 We call a fixed point non-trivial if (αg = 0, αy = 0). 13 We say that saddle-node bifurcation curves with (αg > 0, αy > 0) and (αg > 0, αy = 0) or (αg = 0, αy > 0) are related, if the equilibria that are involved can be related through a transcritical bifurcation in the αg = α nt g = 0 or αy = α nt y = 0 subspace. Figure 20a shows the bifurcations of the sets of fixed points which were indicated as pairs 8 and 10 in the previous section. Here, we see the transcritical and saddlenode bifurcations that were found in the previous section. In addition, there is an extra transcritical bifurcation at the solid red line, where a non-trivial fixed point branches of. The solid red line converges to the branching point at x = 5.5 that was found in the Veneziano limit. Furthermore, we see a transcritical bifurcation at non-trivial value of α y (dash-dotted red line), where a non-trivial fixed point splits of. This transcritical bifurcation converges to the branching point found at x = 11.6 in the Veneziano limit. The two non-trivial fixed points disappear through a saddle-node bifurcation at the solid blue line, which is related to the blue dashed line representing the saddle-node bifurcation that was found in the α y = 0 space. The solid blue line converges to limit point that was found at x = 5.5 in the Veneziano limit. Figure 20b shows the bifurcations of the sets of fixed points which were indicated by the pairs 3 and 4 in the previous section. Here, we see the transcritical and saddlenode bifurcations that were found in the previous section. In addition, there is an extra transcritical bifurcation at the solid red line, where a non-trivial fixed point branches of. The solid red line converges to the branching point at x = 5.5 that was found in the Veneziano limit. Furthermore, we see a trancritical bifurcation at non-trivial value of α y (dash-dotted red line), where a non-trivial fixed point splits of. This transcritical bifurcation converges to the branching point found at x = 11.6 in the Veneziano limit. The two non-trivial fixed points disappear through a saddle-node bifurcation at the solid blue line, which is related to the blue dashed line representing the saddle-node bifurcation that was found in the α y = 0 space. The solid blue line converges to limit point that was found at x = 5.5 in the Veneziano limit. Figure 20c shows the bifurcations of the sets of fixed points which were indicated by the pairs 6 and 9 in the previous section. Here, we see the transcritical and saddlenode bifurcations that were found in the previous section. In addition there is an extra transcritical bifurcation at the solid red line, where a non-trivial fixed point branches of. The solid red line converges to the branching point at x = 7.4 that was found in the Veneziano limit. Furthermore, we see a transcritical bifurcation at non-trivial value of α g (dash-dotted black line), where a non-trivial fixed point splits of. This transcritical bifurcation diverges in the Veneziano limit. The two non-trivial fixed points disappear through a saddle-node bifurcation at the solid blue line, which is related to the blue dashed line representing the saddle-node bifurcation that was found in the α y = 0 space. The solid blue line intersects the black dash-dotted line at SNT6. Therefore only in a small interval N c ∈ (1.2, 1.5) the non-trivial fixed points disappear through a saddle-node bifurcation, while at larger N c the non-trivial fixed points that branch of from the solid red line and the black dash-dotted line diverge for small values of x. In addition, for these fixed points we find a saddle-node bifurcation curve not only in the α y = 0 space, as was the case in the previous figures, but also in the α g = 0 space (dash-dotted blue line). Figure 20d shows a few more transcritical bifurcation curves lying in the α nt g = 0 and/or α nt y = 0 subspace. Non-trivial fixed points that are related to these transcritical bifurcation curves, all diverge in one or multiple couplings for small x instead of disappearing through JHEP07(2019)075   Table 6. Eigenvalues and critical eigenvectors crossing marginality of the saddle-node bifurcation with non-trivial α g , α y at N c = 3.

JHEP07(2019)075
Bifurcation that this complex scaling behavior involves operators, which are linear combinations of the effective four-fermi interactions. Furthermore, we find multiple bifurcation points in figure 20, which are reported in table 7. The first cusp bifurcation and the points SNT5, IP1 and IP2 were already found in the previous model. In addition, we find an extra cusp point in the saddle-node bifurcation curve in the α g = 0 space. This curve has another cusp bifurcation at IP1, making IP1 even more degenerate. The points TCT1 up to TCT6 all describe the intersection of two transcritical bifurcation curves, and can be described by the normal form ẋ = αx + x 2 , y = βy + y 2 . (3.7) Finally, IP3 describes a point where 2 saddle-node bifurcation curves and 2 branches of a transcritical bifurcation intersect, and can therefore be described as a saddle-nodetranscritical bifurcation, where an additional saddle-node bifurcation intersects. Four fixed points are involved in this bifurcation.

JHEP07(2019)075
(a) The probable scenario for the conformal window with completely stable Banks-Zaks fixed point.

Discussion of the results
In this section, we have seen how the addition of a scalar field alters the bifurcations found in the previous section. In particular, we would like to know whether we find similar behavior in the (α g > 0, α y > 0) space as we found in the (α g > 0, α y = 0) space. We find the transcritical bifurcation at x = 5.5, which marks the upper edge of the conformal window. This transcritical bifurcation now indicates the intersection of three fixed points: one on the (α g = 0, α y = 0) subspace, one on the (α g > 0, α y = 0) subspace, and one on the (α g > 0, α y > 0). Furthermore, at this line there are two different transcritical bifurcation curves for different values of (g S , g V , g V 1 , g V 2 ) as was also found in the previous model. One of these bifurcations is shown in figure 20a and the other in figure 20b.
We see that both saddle-node bifurcations that could indicate the lower edge of the conformal window in the model with α y = 0 have a related saddle-node bifurcation in the (α g > 0, α y > 0) space, indicating that the phase transition at the lower edge of the conformal window is similar to the previous section and involves the annihilation of fixed points, followed by a walking behavior of the coupling constant, that changes continuously to a running behavior. However, the non-trivial saddle-node bifurcations are found at a larger value of x than the ones in the α y = 0 subspace, and are triggered by the α g interaction instead of the g S interaction. The resulting estimates for the conformal window are shown in figure 21. This curve is very close but slightly below x = 5.5 at large values of N c and slightly decreases towards smaller N c . Therefore, the conformal window is much smaller in the theory that contains a scalar field in the adjoint representation.

JHEP07(2019)075
Furthermore, we find that the fixed points, which were labeled by 5 and 7 have no related fixed points with non-trivial value of α y . Therefore, we do not have the spiraling behavior related to those fixed points as we had in the previous section. For the saddle-node bifurcation of the pair labeled by 6 and 9 there is only a related bifurcation with non-trivial and positive α g and α y in a small interval N c ∈ (1.2, 1.5), which makes it irrelevant for the physical model, since we require integer values of N c . Therefore, we find that the non-trivial fixed points of the sets 1, 2, 6 and 9 all diverge 14 for small x.

Conclusions and outlook
We have evaluated an effective model for QCD 4 , and studied the fixed points and their bifurcations in and close to the conformal window. We have found a fixed point merger (saddle-node bifurcation) in the Veneziano limit at N f /N c = 4.049 indicating the lower edge of the conformal window. This prediction is similar to the ones by [3,15], where the same model was analyzed. Furthermore, we have continued this bifurcation to the small N c regime, and found that this curve is more or less constant up to (N c = 3.00, N f = 12.17), while For N c ∈ {1, 2} the window is much smaller. On the other hand, there exists another fixed point merger (corresponding to the same merger in the Veneziano limit), which is more or less constant up to (N c = 3.00, N f = 11.26). For N c = 1, these fixed points diverge instead of disappearing through a saddle-node bifurcation, while N c = 2 is an intermediate value where they merge at large coupling constants. Both mergers could in principle indicate the lower edge of the conformal window. Although the second is destabilized by the g V 1 and g V 2 interactions. In both cases the transition at the lower edge of the conformal window is an infinite order phase transition, generated by the effective scalar four-fermi interaction for N c > 2. Furthermore, we have found a few extra sets of fixed points, that either disappear through saddle-node bifurcations or diverge. We stress that the saddle-node bifurcation is stable against small perturbations of our model, while the divergence scenario is not.
The existence of multiple fixed points and multiple saddle-node bifurcations in this model, opens up many new questions about the lower edge of the conformal window. The infinite order perturbative beta function may contain more fixed points apart from the well known Banks Zaks fixed point. This could result in a one or more transitions within the conformal window, where pairs of fixed points disappear through mergers. The lower edge is then reached when the last fixed point merger takes place. Although the model is stable against small deformation, it would be interesting to extend the model by including running coupling constants for higher dimensional effective interactions, and see which of the fixed points and saddle-node bifurcations remain to exist.
In this paper, we already extended the effective QCD model with an scalar field coupling through a Yukawa interaction to the fermions. This interaction turns out to destabilize the model. In this model we found that the 2 candidates for the conformal window both exist in this model as well. Although the lower edge of the conformal window is JHEP07(2019)075 located at larger values of N f /N c for such a model, as was already predicted in [34]. The other fixed points, that were found in the model without Yukawa interactions, are pushed into a regime α g < 0 and/or α y < 0.
In addition, we have found that the model allows for complex scaling dimensions of the effective interactions. These complex scaling dimension contain both a scaling and a rotation of the associated operator. The rotation indicates that the (ir)relevant operators associated with the complex eigenvalue rotates in theory space, while the orbits move away from/towards the fixed point. Such behavior was found earlier in for example [36,37]. For one pair of fixed points in the effective QCD model these complex scaling dimensions induced a clearly visible spiraling behavior in the RG-flow. This pair disappeared in the model extended with a scalar field and Yukawa coupling. It would be interesting to study these kind of exotic orbits in more detail, and maybe relate them to other exotic RGflows that have been found, such as the bouncing solutions recently found in holographic studies [38].
Finally, we have shown that the use of bifurcation analysis in the study of RG-flows can be very useful tool in mapping out the fixed points and their transition in complicated model with multiple variables and parameters as was suggested in [3]. It would be interesting to use numerical tools such as the MATLAB package MatCont [1,2] in other studies of renormalization group flows. In particular for non-perturbative beta functions this can lead to new insights. An alternative to our numerical approach would be to use Conley index theory and describe the exit sets for RG-flows as was suggested in [3]. This would then lead to a less detailed picture, from which already part of the relevant information can be deduced.
Furthermore, it could be interesting to use general statements about relations between symmetries in the action and symmetries of the beta functions together with bifurcation theory to determine whether certain bifurcations can be expected in specific physical models. This could help to determine whether a specific model can have fixed points based on the symmetries of the model.

JHEP07(2019)075
typically T = N 0 or T = N for discrete dynamical systems, and T = R + 0 or T = R for continuous dynamical systems. The collection of all maps {Φ t } t∈T is called the flow of the dynamical system and should satisfy Therefore, the flow forms a group, if Φ is invertible, and a semi-group otherwise.
For a point x ∈ X, one can define the orbit O Φ (x) := {Φ t (x)|t ∈ T )}. Furthermore, one can define the notion of an invariant set of a dynamical system.
Definition. Given a dynamical system (T, X, Φ t ), • a set that is both forward and backward invariant is called an invariant set of the dynamical system.
If the map Φ t is non-invertible and T is only defined on positive numbers, every forward invariant set is called an invariant set of the dynamical system. It is assumed above that Φ t (x) is defined for all (t, x) ∈ T × X. If it is not the case, i.e., given x ∈ X, Φ t (x) is only defined for a subset of T that includes t = 0, the dynamical system is called local, and all definitions should be modified accordingly.
A natural way to generate a dynamical system is by a set of ordinary differential equations (continuous-time dynamical system) or a set of difference equations (discrete-time dynamical system). In this article, we consider dynamical systems generated by a set of differential equations, since those can be related to the beta functions of the renormalization group. Given the space X = R n , one can consider the system of autonomous ordinary differential equationsẋ = f (x), (A.1) where x ∈ R n and f : R n → R n , which is assumed to be (sufficiently) smooth. If one defines the map Φ t : R n → R n by Φ t (x 0 ) = x(t, x 0 ), then the tuple (R, R n , Φ) is a (local) continuous-time dynamical system. The solutions for specified initial conditions, x(0, x 0 ) = x 0 , define the orbits and the visualization of the flow is called the phase portrait. In the study of dynamical systems equilibria 15 and periodic, homoclinic, and heteroclinic orbits are of particular interest.
Definition. Given a continuous-time dynamical system (R, R n , Φ t ), then a point x ∈ R n is called an equilibrium if Φ t (x) = x ∀t ∈ R.

JHEP07(2019)075
Definition. Given a continuous-time dynamical system (R, R n , Φ t ), a point x ∈ R n , and an orbit O Φ (x), then the orbit is called • a periodic orbit if x is not an equilibrium and Φ T (y) = y ∀y ∈ O Φ (x) for some T > 0; minimal such T is called the period; where y is an equilibrium; The topology of the phase portrait is characterized by its invariant sets and their stability. An invariant set is asymptotically stable, if all orbits starting in a neighborhood of the invariant set stay in this neighborhood and converge towards the invariant set if t → ∞ and unstable otherwise. For generic equilibria the stability is determined by the eigenvalues of Jacobian matrix at the equilibrium. At the equilibrium, one can define the notions of a stable, unstable and critical eigenspace.
An equilibrium is called hyperbolic if dim(T c ) = 0, stable if dim(T s ) = n and unstable if dim(T u ) ≥ 1. Furthermore, an equilibrium with dim(T s ) = n or dim(T u ) = n is a node and an equilibrium with dim(T s ) ≥ 1 and dim(T u ) ≥ 1 is a saddle. The classification of non-hyperbolic equilibria not satisfying these conditions is more involved, and will not be discussed further in this article. For a generic equilibrium in a continuous-time dynamical system, one can also define the notions of a stable and unstable manifold of the equilibrium.
Definition. Given a continuous-time dynamical system (R, X, Φ t ), and an equilibrium point x ∈ X of Φ, then • the stable manifold of x is defined as W s (Φ, x) = {y ∈ X : lim t→∞ Φ t (y) = x}, and • the unstable manifold as W u (Φ, x) = {y ∈ X : lim t→−∞ Φ t (y) = x}.

JHEP07(2019)075
For a dynamical system defined by (A.1), the stable manifold is at the equilibrium tangent to the stable eigenspace, while the unstable manifold is at the equilibrium tangent to the unstable eigenspace. For nonhyperbolic equilibria, one can introduce the concept of a center manifold W c loc (Φ, x), which is n c -dimensional and tangent to the critical eigenspace at the equilibrium. The existence of this locally defined manifold is implied by the Center Manifold Theorem.

A.2 Bifurcation analysis
A continuous-time dynamical system defined by a set of ordinary differential equations can be analyzed by solving the ODE's numerically. However, if the differential equations also depend on a set of parameters, f α (x), α ∈ R m , one is often not interested in the exact solutions for specific values of α, but in the topological properties of the phase portrait such as the number of equilibria, limit cycles and their stability or the existence of homoand heteroclinic orbits and chaotic attractors. In general, changes of the phase portrait, when α is varied, occur smoothly except for some specific values of α, where the topology of the phase portrait changes. Changes in the topology of the phase portrait are called bifurcations. Topology changes often occur near or within invariant sets. Bifurcations can be divided in 2 subclasses: • local bifurcations: bifurcations that can be detected by considering small but nonshrinking neighborhoods of equilibria or limit cycles.
• global bifurcations: bifurcations that cannot be detected in this way. Topology changes typically involve multiple invariant sets.
In addition, bifurcations may be called subcritical or supercritical, depending on the stability of the equilibria or limit cycles that appear or disappear at the bifurcation. Furthermore, bifurcations can be characterized by their codimension. The codimension of a bifurcation in a generic system is given by the number of independent conditions that have to be satisfied for a bifurcation to happen, and for local bifurcations is often related to the number of critical eigenvalues of the Jacobian matrix at the equilibrium. Consequently, the codimension of a bifurcation is always smaller or equal to the dimension of the parameter space. In general, many types of bifurcations can be found. For an overview of possible bifurcations and their conditions, we refer to [40]. In this appendix, we discuss bifurcations and related terminology that is relevant for this article. In doing so, we closely follow [40]. Bifurcations in low dimensional systems and with low codimension can often be found using analytical techniques, but for higher dimensional systems and in particular for higher codimension of the bifurcation one often has to rely on numerical techniques. Several computer programs have been developed for bifurcation analysis. In this article we make use of the Matlab package Matcont [1,2]. Bifurcation software relies on the conditions that have to be satisfied in order for a bifurcation to happen.
In different dynamical systems one can expect different bifurcations, which are not immediately defined by the same conditions. However, one would like to have a way of JHEP07(2019)075 categorizing bifurcation in various systems. For this, we need the notion of topological equivalence. Consider two continuous-time dynamical systems: Since many bifurcations are locally defined, one is often only interested in local topological equivalence near specific points.
Definition. Dynamical systems (A.2) and (A.3) are locally topologically equivalent near the origin, if there exists a map (x, α) → (h α (x), p(α)), defined in a small neighborhood of (x, α) = (0, 0) ∈ R n × R m such that • p : R m → R m is a homeomorphism defined in a small neighbourhood of α = 0, p(0) = 0; • h α : R n → R n is a parameter-dependent homeomorphism defined in a small neighborhood U of x = 0, h 0 (0) = 0, and mapping orbits of the system (A.2) in U onto orbits of the system (A.3) in h α (U ), preserving the direction of time.
By coordinate translations one can easily generalize this definition to local topological equivalence near arbitrary points. Having these definitions, one can categorize bifurcations by topological equivalence to normal forms. Consider a polynomial continuous-time dynamical systemż = g(z, β; σ), z ∈ R n , β ∈ R k , σ ∈ R l , (A. 4) which has an equilibrium z = 0 at β = 0, satisfying k bifurcation conditions. Here, σ is a vector of the coefficients of the polynomial g(z, β; σ).
Definition. System (A.4) is a topological normal form for the bifurcation if any generic system (A.2) with equilibrium x 0 satisfying the bifurcation conditions α 0 is locally topologically equivalent near (x 0 , α 0 ) to (A.4) for some values of the coefficients σ i .
We notice that the system (A.2) can have higher dimension than is needed for the bifurcation to occur: x ∈ R n , while the normal form has z c ∈ R nc with n c < n. In such a case, there is a continuation of the critical center manifold. where g(z c , β; σ) is a normal form on the center manifold and z s ∈ R ns , z u ∈ R nu correspond to the stable and unstable eigenspaces. We now give a short discussion of local bifurcations which are encountered throughout the article. We begin with bifurcations in scalar one-parameter ODEṡ In the n-dimensional case, each normal form below is the normal form (A.5) on the onedimensional center manifold.
Saddle-Node bifurcation. This bifurcation is a codimension one bifurcation at which two equilibria collide and disappear, and is also known as fold or limit point bifurcation.
Transcritical bifurcation. This bifurcation is a codimension one bifurcation in the class of systems that always have a trivial equilibrium, and is also a branching point: at the bifurcation, two equilibria collide and exchange their stability. The bifurcation occurs at (x,ᾱ), if the following conditions hold for the system (A. • f (x,ᾱ) = 0, The system is then locally topologically equivalent near (x,ᾱ) to the normal forṁ y = βy ± y 2 , y, β ∈ R (A. 10) near (y, β) = (0, 0) (see figure 23). Notice that we can transform the normal form of the tanscritical bifurcation to the normal form of the saddle-node bifurcation by the noninvertible coordinate transformation This is called a nonversal unfolding of the saddle-node bifurcation.
Pitchfork bifurcation. This bifurcation is a codimension one bifurcation for systems with reflectional symmetry, and is also a branching point: at the bifurcation one equilibrium splits into three equilibria. The bifurcation occurs at (x,ᾱ), if the following conditions hold for the system (A.8): • f xα (x,ᾱ) = 0.
The system is then locally topologically equivalent near (x,ᾱ) to the normal forṁ y = βy ± y 3 , y, β ∈ R (A.11) near (y, β) = (0, 0) (see figure 24). The − sign gives the supercritical case where a stable equilibrium splits into two stable equilibria and an unstable equilibrium, and the + sign is the subcritical case, where an unstable equilibrium splits into two unstable equilibria and a stable equilibrium. We continue with two bifurcations in scalar two-parameter ODEṡ In the n-dimensional case, each normal form below is the normal form (A.5) on the corresponding one-dimensional center manifold.
Cusp bifurcation. This bifurcation is a codimension two bifurcation in generic systems, where two branches of saddle-node bifurcation meet tangentially. Nearby, the system can have three equilibria which collide and disappear at the saddle-node bifurcations. The bifurcation occurs at (x,ᾱ), if the following conditions hold for the system (A.12): • f (x,ᾱ) = 0, where c is the cubic coefficient of the system (A.12). The system is then locally topologically equivalent near (x,ᾱ) to the normal forṁ y = β 1 + β 2 y ± y 3 , y, β , β 2 ∈ R (A. 13) near (y, β) = (0, 0) (see figure 25). Saddle-Node-transcritical bifurcation. This bifurcation is a codimension two bifurcation in the class of systems that always have a trivial equilibrium, where two branches of a transcritical bifurcation curve meet on a saddle-node bifurcation curve. There exist two types of this bifurcation. One can be seen as a nonversal unfolding of a Cusp bifurcation, while the other can be obtained as an unfolding of a degenerate Bogdanov-Takens bifurcation. Here, we discuss only the first case. Nearby, the system can have three equilibria: one pair has a transcritical bifurcation at one of the branches and another pair has a transcritical bifurcation at another branch. The last possible pair of equilibria collides and disappears at the saddle-node bifurcation curve. If the bifurcation occurs at (x,ᾱ), then the system (A.12) is locally topologically equivalent near (x,ᾱ) to the normal form given bẏ y = β 1 y + β 2 y 2 ± y 3 , y, β 1 , β 2 ∈ R (A.14) near (y, β) = (0, 0) (see figure 26). We then find a transcritical bifurcation along the line β 1 = 0, and a saddle-node bifurcation along β 1 = ± β 2 2 4 . The normal form can be transformed to the normal form of the cusp bifurcation through the transformation Finally, we discuss for completeness the last codimension one local bifurcation in generic ODEs (A.2), i.e. Andronov-Hopf bifurcation. Note that it does not occur in our RG-flows but might appear in other models. For this bifurcation, the center manifold is two-dimensional, so it is sufficient to consider (A.2) with n = 2 and m = 1.
Andronov-Hopf bifurcation. This bifurcation is a codimension one bifurcation, and implies the birth of a periodic orbit (limit cycle). The bifurcation occurs at an equilibrium (x,ᾱ), if a pair of two conjugate complex eigenvalues of the Jacobian matrix crosses the JHEP07(2019)075  Figure 26. Saddle-node-transcritical bifurcation in the systemẏ = β 1 y + β 2 y 2 + y 3 .

B.1 Feynman rules for QCD 4 with a four-fermi interaction
The theory as discussed in section 2 has three propagators: 16 where µ, ν label the space-time indices, A, B label the colors in the adjoint representation of SU(N c ), a, b label the color in the fundamental representation of SU(N c ) and {i, j} label the flavor in the fundamental representation of SU(N f ). Furthermore, we take the Landau gauge in which the gauge parameter ξ = 0.
The fermion-gluon interaction (gauge vertex) is given by: 16 All Feynman diagrams are generated with the Tikz-Feynman package [41].

JHEP07(2019)075
Furthermore, there are four four-fermi interactions: Notice that the vertices have been split up in two parts. This is done to make clear along which lines the color charge is conserved, and will make the calculations in appendix B.3 more tractable, but has no further physical meaning.

B.2 Beta function for the gauge vertex
In this appendix, we reproduce the method used in [15] to calculate the beta function for the gauge vertex. We're only interested in gauge invariant contributions to the beta function, but the Wetterich equation does not necessarily respect gauge invariance. Therefore, we do not use the exact renormalization approach. Instead we start from the well known 2-loop JHEP07(2019)075 perturbative beta function for QCD 4 with N c colors and N f massless flavors [5]: Next, we would like to include effective interactions induced by the four-fermi couplings. These can be added perturbatively to the beta functions. The beta function is then given by where the last term is represents the perturbative corrections induced by the effective couplings. These perturbations should be chosen such that they're gauge invariant and such that they represent perturbations, that are included in lim n→∞ β [n] g , but not in β [2] g . There is one one-loop correction induced by the effective four-fermi interaction, which is represented by = However, this correction turns out to be gauge dependent [31], and therefore we discard it. A possible two-loop correction is given by = For the effective four-fermi interaction one could take O i , i ∈ {S, V, V 1 , V 2 }. However, since the interactions should represent perturbations, which are included in lim n→∞ β [n] g the chirality should be the same in the whole diagram. This excludes the interactions O S and O V 1 . Furthermore, due to the tracelessness of the generators of SU(N ) the contribution due to O V 2 vanishes. Therefore, only O V contributes, and its contribution is in the sharp cutoff limit given by yielding a contribution to the beta function of the form

B.3 Beta functions for the four-fermi interactions
We can find the beta functions of the model using the equation [42] d dt Furthermore, For this Lagrangian the anomalous dimension is η = 0 [15]. The last term can be evaluated using the Wetterich equation [43] dΓ The right hand side is equal to the expansion in terms of all 1-particle-irreducible diagrams, so we only have to find all the 1-particle-irreducible diagrams proportional to the four point interactions. We'll evaluate the diagrams in the sharp cutoff limit, where δR Λ = Λδ(|p| − Λ)δΛ. Furthermore, we use spherical coordinates in which where the Ω represents the spherical part. If there is no angular dependence this part is the surface area of a 3-dimensional sphere, which is 2π 2 . In the calculations of the diagrams JHEP07(2019)075 we'll make use of the following integrals, which follow from symmetry arguments dp p µ p ν p 2n = In addition, we make extensive use of the following identities Tr(γ µ γ ν ) = 4η µν Tr(γ 5 ) = Tr(γ µ γ ν γ 5 ) = 0, ε αµνρ ε βµνρ = −6δ β α .
Using the Feynman rules from appendix B.1, we find the following contributing 1PI diagrams: = + +

JHEP07(2019)075
There are 29 different diagrams of the first type, 12 of the second type and 4 of the third. Five diagrams of the first type contribute to the beta function of G S . One of those is proportional to g 2 S : which is evaluated as One diagram is proportional to g S g V : which is evaluated as 17 N f G S G V 2π 2 Λ 2 O S . Two diagrams are proportional to g S g V 1 : 17 Due to lengthiness of the calculations, we have suppressed further explicit calculations of the diagrams in this article.

JHEP07(2019)075
One diagram is proportional to g S g V 2 : which is evaluated as G S G V 2 2π 2 Λ 2 O S . Five diagrams of the first type contribute to the beta function of G V . One of those is proportional to g 2 S :

JHEP07(2019)075
Nine diagrams of the first type contribute to the beta function of G V 1 . One of those is proportional to g 2 S :

JHEP07(2019)075
Two diagrams are proportional to g V g V 1 : which are evaluated as Two diagrams are proportional g V 1 g V 2 : which are evaluated as Ten diagrams of the first type contribute to the beta function of G V 2 . Two of those are proportional to g 2 V : which are evaluated as One diagram is proportional to g 2 V 1 :

JHEP07(2019)075
which is evaluated as Four diagrams are proportional to g 2 V 2 : which are evaluated as One diagram is proportional g S g V 1 : which is evaluated as which are evaluated as

JHEP07(2019)075
Three diagrams of the second type are proportional to g S α g and contribute to the beta functions of G S and G V 1 : The second and the third both evaluate to 0 and the first is evaluated as Three diagrams of the second type are proportional to g V 1 α g and contribute to the beta functions of G S and G V 1 : The first and second both evaluate to 0 and the third evaluates to Three diagrams of the second type are proportional to g V α g and contribute to the beta JHEP07(2019)075 functions of G V and G V 2 : The first and third evaluate to 0 and the second is evaluated as Three diagrams of the second type are proportional to g V 2 α g and contribute to the beta functions of G V and G V 2 : The first and third evaluate to 0, while the second is evaluated as Furthermore, the third type of correction is proportional to α 2 g and consists of 2 diagrams contributing to the beta functions of G S and G V 1 :  32 1 and two contributing to the beta functions of G V and G V 2 : which are evaluated as

C Beta functions for QCD 4 with a scalar field
When the scalar terms (3.2) are added to the model (2.4), the set of beta functions is changed. As a consequence, new Feynman rules are added: -60 -

JHEP07(2019)075
First, we consider the function β αg . Due to the addition of the scalar field terms, the two loop beta function for β αg is changed to [44] β [2] g ≡ Λ The additional terms due to the four-fermion interactions can then be added as was done in appendix B.2. The beta function β αy is given by [45] Λ dα y dΛ = 2α y γ φ + γψ ψ .
The anomalous dimension of the scalar field is given by 2N c α y [44]. Using the four-fermi interactions, one can find a non-perturbative approximation to the anomalous dimension of the mass of the fermion field [34] such that we get three contributing diagrams + + .
The first two diagrams contribute N f α y and −6 N 2 c −1 2Nc α g , and two diagrams can be found that give a contribution to the second diagram: which contribute −2N c g S and 8g V 1 .
Next, we consider the beta functions four the four-fermi interactions, which are a slight variation of the ones found for the model without scalars in appendix B.3. The anomalous dimension of the fermion fields is changed due to the presence of the scalar fields to η = α y [34]. Furthermore, we find additional contributions to represented by three types of diagrams + + .

JHEP07(2019)075
Diagrams represented by the third type all evaluate to 0. For the first type we find four diagrams that each contribute to one of the beta functions of the four-fermi interactions: contributing respectively 4g V α y to β g S , g S α y to β g V , −2g V 2 α y to β g V 1 and −2g V 1 α y to β g V 2 .
Furthermore, for the second type there are 2 contributing diagrams: contributing 2α 2 y to β g V 1 and 2α 2 y to β g V 2 respectively.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.