Disordered fermions, extra dimensions and a solvable Yang-Mills theory

Generalizing disorder couplings of the SYK model by means of SU(N) matrices we formulate a lattice model of fermions in d+1 dimensions. Integration of fermions yields an effective theory of Yang-Mills fields in d dimensions, the latter approaching the standard Yang-Mills theory in the case of heavy fermions and the classical limit of vanishing coupling constant of the theory. Quantum mechanically, the theory is solved using large N approximation of the dual effective theory of Hermitian matrices in d dimensions. The theory is asymptotically free and confines the color. In case of massless fermions the emerging theory is an asymptotic safe QCD theory. We discuss also the relationship of this theory to the SYK model.


Introduction
Yang-Mills theories and Quantum Chromodynamics (QCD) [1][2][3][4] play an important role in our understanding of basic forces of the Universe. They are part of the Standard Model of particle physics. Nonetheless, an analytical solution in four dimensions is missing. Although early lattice simulations of Creutz were illuminating and showed that there is no second order phase transition between strong and weak coupling [6], the possibility of a Gross-Witten and Wadia transition [7,8] diminished the hope of a full solution from the strong coupling regime. While an analytical solution is highly desirable, the high precision lattice simulation of Lüscher and Weisz leaves no doubt that, at low energies, the Yang-Mills theory may be described by an effective string theory [9]. Many years of research in lattice simulations of QCD have proven to be an indispensable tool in understanding the Standard Model at a non-perturbative level. During these years, string theory has made a great contribution in the search for a unifying theory of gravity and particle physics. The AdS/CFT correspondence, put forward first by Maldacena [10], has established an avenue in this direction. There is recent progress made in the field using the generic model of Sachdev-Ye and Kitaev (SYK) [11][12][13], which uses disorder and extra dimensions as model building ingredients.
In this paper, we formulate a solvable model beyond the Standard Model. The basic degrees of freedom of the model are d + 1 dimensional lattice fermions coupled by means of SU(N) random matrices defined on a d dimensional sublattice. Such matrices generalize the randomly distributed couplings of the SYK model. However, such a generalization has a non-trivial effect in the structure of the theory: the model is chosen with the aim of obtaining a lattice gauge theory as an effective theory which remains after fermion integration. There may be an arbitrary number of such models. We have constrained the model with the intention to have a dual description of gauge theories and quantum mechanical models of holography such as SYK model and its generalizations.
Note that one may use complex N × N matrices as disorder fields. In appendix A we formulate a matrix theory with Gaussian disorder couplings. In this approach, SU(N) matrices embedded in general complex matrices are the basis of an emergent Yang-Mills theory.
In our model, we fix the length of the extra dimension by the coupling constant of the theory. Therefore, the d dimensional system evolves along the extra dimension in order to reach a certain value of the coupling constant. Since each value of the latter corresponds to a given length scale, the extra dimension in our theory is the dimension of physical scales.
The results of this paper may be summarized in the following: • the model with heavy fermions yields an effective theory of local Yang-Mills fields, a theory which approaches the standard Yang-Mills theory in the classical limit of vanishing coupling constant; • the model with massless fermions yields an effective theory of QCD with a large number of light fermions; • in the large N limit the theory is dual to a field theory of matrices of order N t , where N t is the number of lattice sites along the extra dimension; • the theory is non-perturbatively solvable in the large N approximation; • the theory with heavy fermions is asymptotically free and confines the color; • the theory with massless fermions is dual to an asymptotic safe QCD theory, i.e. a QCD theory which is ultraviolet complete at a non-zero value of the coupling; • the massless fermion theory is also chirally symmetric, a symmetry which spontaneously broken; • although the heavy fermion theory shares important properties with the standard Yang-Mills theory like asymptotic freedom and color confinement, the corresponding beta functions are different functions of the coupling constant; • the asymptotic safe QCD theory, within the leading approximation of our bosonization approach and the low energy limit, is equivalent to an ideal gas of q = 2 SYK models.
In summary, we show that a local Yang-Mills theory is solvable in the dual formulation in the large N limit. Its solution shares qualitative properties of the standard Yang-Mills theory. Moreover, we are able to link gauge theories with the SYK model of holography.
The paper is organized as follows: in the next section we define the model and discuss the disorder average of the theory. In section 3 we show that the effective description of the theory with heavy and massless fermions. In the first case we obtain a local theory Yang-Mills fields in the weak coupling limit, whereas in the second we get QCD with a large number of light modes. Section 4 deals with the disorder averaged theory which is a field theory of matrices of order N t . In section 5 we solve the matrix field theory in the large N limit and heavy fermions. Section 6 deals with the case of massless fermions in which case an asymptotic safe QCD theory emerges. In section 7 we discuss the relation to the q=2 SYK model. The last section closes the paper with the discussion of results and their phenomenological relevance within and beyond the Standard Model. The paper includes three appendices. In appendix A we show that the model can be generalized in the case of Gaussian disorder, where we get an identical large N solution as with SU(N) disorder within our approximation of bosonization of the fermion theory. In order to make the paper self contained, appendix B serves as a starting point in group integration and one-link integrals. Finally, appendix C applies our bosonization approach to strong coupling QCD.

The model
In this section we define our model using Hamiltonian and field theory formalism. We begin with the description of the Hamilton operator.

Hamilton operator
Let Ψ α c (x), Ψ α c (x) * , c = 1, 2, . . . , N be N Dirac-fermion annihilation and creation operators for each Dirac component α = 1, 2, . . . , 2 d/2 at each site x = (x 1 , x 2 , . . . , x d ) on a regular Euclidean lattice of even dimensions d acting on the Hilbert space H. They obey the anticommutation relations: The lattice is finite and we assume it to be a torus with V number of sites. Let also U µ (x) be a SU(N) matrix at each directed link (x, x +μ) on the lattice, where µ = 1, 2, . . . , d. The Hamiltonian operator of the model is: where a is the lattice spacing, κ a dimensionless coupling constant, γ µ are the usual gammamatrices and m is the bare fermion mass. We have denoted by γ 5 the product of all gammamatrices in order to remember that the theory in four dimensions is the phenomenologically relevant theory. In the above expression we have suppressed Dirac indices for clarity. In the following the lattice spacing is set to unity. Note also that we have chosen naive fermions on the lattice. Any local discretization would do the job. For example, the Kogut-Susskind version [4]: 1 is simpler since fermion operators carry no Dirac indices, i.e. they obey the anticommutation relations: Here,γ 5 is the lattice site parity operator taking ±1 values on even/odd lattice sites, i.e. γ 5 (x) = (−1) x 1 +···+x d and η 1 (x) = 1, η µ (x) = (−1) x 1 +···+x µ−1 , µ = 2, . . . , d. Both,γ 5 and η µ are diagonal matrices on the lattice. If we define hopping matrices: then the following commutation relations hold: For the clarity of notations we introduce the Hermitian fermion matrix: This way, our model may be written in the form: where h(x, y) cc ′ are the matrix elements of the V N × V N fermionic matrix h. Note that the Hamilton operator (2.3) defines a dynamical system in d dimensions, where d is the number of dimensions of the world we live in, i.e. d = 4. The system evolves along an extra dimension, which for the moment may be thought as a fictitious time dimension. We take this extra time to be imaginary, which means that we are studying the system at a finite temperature with partition function: where the trace is taken in the Hilbert space and N t is the length of of the extra dimension in units of the lattice spacing a. Using the Hilbert space trace rules we get: (2.10) In this paper we are interested in the large N t limit of the theory. In the next subsection we compute the effective gauge action of the theory in this limit.

Effective gauge action
Since h is traceless and using the identity det A = e Tr log A we can write the right hand side of (2.10) in the form: where the trace is taken in the tensor product space of the lattice sites and the SU(N) group. The matrix in the last expression is an even function of h, therefore we get: (2.12) In the large N t limit the second exponential may be neglected if mN t is not smaller than O(1). 2 With this assumption, the effective action of the theory is: Eventually, we would like to integrate gauge fields and the right hand side is not convenient for this purpose. In the next subsection we formulate the theory as a d+1 field theory suitable to deal with the integration of the SU(N) field.

Fermion action
Another way to formulate the model is by introducing Grassmann valued fermion fields Θ(x, t) andΘ(x, t), where t = 1, 2, . . . , N t label the lattice sites along the extra dimension. N t assumes integer values and the Grassmann field satisfies antiperiodic boundary conditions along the extra direction. The action of the theory is: where∂ t is a lattice derivative, for example: In continuum limit this formulation of the model is equivalent to the Hamilton operator formulation (2.3).

Gauge invariance
The action (2.14) defines a lattice theory of a massive fermion in d + 1 dimensions in a fixed background of SU(N) random matrices, which are defined on the links of the d dimensional sublattice replicated along the extra dimension. As it stands, it allows gauge variant states to propagate in the bulk. In order to ensure a full gauge invariance one can introduce an extra SU(N) field along the extra dimension and different random SU(N) matrices at each time slice of the theory, i.e. U µ (x, t), µ = 1, 2, . . . , d, t = 1, 2, . . . , N t . Fixing the axial gauge, the gauge invariant action within the restricted class of gauge transformations of the axial gauge is given by: where the extra dimension SU(N) field resides at the boundary. Note also that there are N t different fermion matrices at each time slice of the lattice: where we have suppressed space and color indices for clarity (recall equation (2.5)). This model differs from the model (2.14) in two ways: it is gauge invariant along the extra dimension and gauge fields along this dimension are t dependent. The difference between these two models should be clarified in a separate study. In the following we focus with the time independent gauge field model (2.14).

Action in canonical form
In order to bring the action in the form of a d + 1 Dirac theory we absorbγ 5 in a new set of fields: In terms of these fields the action of equation (2.14) takes the form: We use this action in the following. Note that we have assumed an even number of dimensions d in order to be compatible with d = 4 dimensions of the Standard Model in the case of infinite N t limit. However, for the sake of generality one may allow d to be an odd integer too. In this case the Dirac field has (d + 1)/2 components and the field ψ(x, t) may not necessary be an irreducible spinor representation. Moreover, Kogut-Susskind fermions may be formulated in any dimensions.

Disorder average
Theories with disorder such as SYK use the replica trick and compute disorder average by integrating the Boltzman kernel of the replicated theory with respect to the probability distribution of the random couplings. In the limit of large N it is customary to assume a saddle point solution of decoupled replicas. Gross and Rosenhaus as well as Kitaev and Suh discuss this issue in the appendix of a recent publication on the SYK model [13,14]. The net result of the diagonal replica assumption is that the disorder averaged theory is, in the large N limit, the disorder averaged partition function of the fermion theory. This result motivates us to promote SU(N) matrices as a randomly distributed field, in which case the partition function of the theory is: where dU µ (x) is the Haar measure of the SU(N) group integration and dψ(x, t) a (as well as dψ(x, t) a ) is the Grassmann integration measure. This definition allows one to identify two effective theories depending which fields are integrated first. Integration of fermion fields gives a SU(N) gauge field theory in d dimensions, while the integration of SU(N) fields gives an interacting fermion theory in d + 1 dimensions. Thus, we have in principle a dual description of the same theory from the start.

Meaning of the extra dimension
We have defined a quantum mechanical model in d dimensions, keeping in mind that d = 4 is the physically interesting case, where one of dimensions is the time dimension. The length of the extra dimension N t is related to the value of the coupling constant κ, as it is shown in the next section. Therefore, the d dimensional system evolves along the extra dimension in order to reach a certain value of κ. Since each value of κ corresponds to a given length scale, the evolution of the system along this dimension is the evolution to reach that scale. This way, the extra dimension is also the dimension of physical scales.
In the next section we show that the d dimensional effective theories derived from the model are of the Yang-Mills and QCD type.

An effective theory of local Yang-Mills fields
In this section we are interested in computing an effective theory of Yang-Mills fields. In the next subsection we proceed with the case of massive fermions.

Massive fermions
The largest fermion mass on a lattice is proportional to m/a and we fix it to be exactly 1/a. Then, form (2.7) we have: This way, using (2.13), the effective action of the theory is: Expanding in powers of κ 2 the right hand side we get: where c o is real, c 1 = 1/4 and we have used the fact that the product of η matrices around the plaquette equals −1. If we want our theory to be the standard Yang-Mills theory in the classical limit of vanishing coupling constant, then we have to make sure that the plaquette term is the Wilson discretization of the standard Yang-Mills theory. Therefore, scaling the length of the extra dimension according to the relation: we get an effective theory which is described by the action: When no other loops are present other than the plaquette, the theory is the Wilson discretization of the standard Yang-Mills theory. Therefore, the above theory is a theory of Yang-Mills fields enlarged by larger Wilson loops. The theory is local if the series expansion in terms of Wilson loops converges. An upper bound of the convergence radius of the series may be found using the following standard argument: since at each lattice site there are 2d − 1 choices to construct a loop without moving backwards, the number of Wilson loops of length n is bounded by (2d − 1) n . This way, at each order n in the κ expansion, the number of Wilson loops of length n does not grow faster than (2d − 1) n . Hence, the series converges for κ ≤ 1/(2d − 1) and the theory is local for vanishing κ.
Classically, as κ goes to zero, the theory is dominated by the plaquette term while larger Wilson loops may be included in a perturbation expansion. The perturbative approach in the classical theory is important if we want to maintain the correspondence principle, in which case, our theory includes the standard Yang-Mills theory as a leading order approximation of the expansion. Quantum mechanically, we may expect that the theory is asymptotically free, since larger Wilson loops are perturbations of the plaquette result. This is indeed the case as shown in section 5. However, the beta function of the theory vanishes linearly with the coupling constant, which shows that the approach to the continuum limit is different from what we expect in the standard Yang-Mills theory.
Note that we have related the length of the extra dimension to the coupling constant of the gauge theory. Given N t is discrete, the relation (3.4) tells us that κ is discrete too. However, since we are interested in the large N t limit of the theory, we assume κ values to be real.

Massless fermions
In case of massless fermions the fermion matrix is κh o . Nonetheless, we introduce a small mass m of the order O(1/N t ). The square of the fermion matrix in this case is: In order to discover an effective local theory we start with the action of the theory (2.14). By means of Fourier transformed fieldsΘ(x, ω k ) andΘ(x, ω k ) the action may be be written in the form:Ĩ where the choice: respects the boundary condition along the extra dimension. Integrating Grassmann fields the partition function may be written as a product of the following determinants: Assuming N t is even, we get: 3 The effective action of the theory: is thus a theory N t /2 Kogut-Susskind fermions (modulo factor four) with masses: Note that the theory has a large number of heavy flavors with masses of the order 1/κ. For such masses one can write an effective action similar to the convergent expansion (3.5). These contribute to the Yang-Mills part of the effective action. The rest of masses m k contribute to the fermionic part of the action. This way, we end up with a local quantum field theory, which is QCD with a large number of flavors. Since a large number of frequencies scale like 1/N t a large number of fermions are light with masses of the order κ 4 (for m ∼ 1/N t ). We expect asymptotic freedom, which is a property of QCD with a limited number of massless flavors, to be lost in this case. In section 6 we show that the theory is asymptotic safe at a fixed value of κ which corresponds to a large value of N t .

Synthesis
In this section we have shown that, when fermions are integrated out, the large N t limit of our model gives effective Yang-Mills and QCD theories. As shown in the rest of the paper, the dual representation of the model is solvable in the large N limit. This property gives us an analytical tool to analyze these theories in this limit. In the next section we deal with the effective theory in terms of dual degrees of freedom.

A Hermitian matrix field theory
The aim of this section is to formulate the theory in terms of dual fields. We first integrate gauge fields obtaining a pure fermion theory with the action S defined by the equation: Then we bosonize fermions in terms of matrix fields. We begin with the gauge field integration.

Integration of gauge fields
The partition function (4.1) may be written as a product of one-link integrals: Group integration rules and one link integrals have been solved long ago by different authors, see for example [16][17][18][19] and references therein. In order to make the paper self contained we have derived one-link integrals in appendix B. The result is: where the solution F (.) depends on a set of boundary conditions (see (B.12)). Since we are concerned with the behavior of the theory in the large N t (i.e. small κ) limit, we require the leading order expansion of F (.) in the form: where, given a lattice site x and direction µ, Λ is a matrix with matrix elements: The full result is (see (B.17)): In this paper we use the leading order approximation (4.4), i.e. the action that we use in the following is: (4.7) In order to estimate the effect of higher order terms we compute the expectation value itself is expected to follow the behavior of Λ(t 1 , t 2 ) in this limit. The expectation value is computed using equation (4.27): where G(x, t, t ′ ) is the two point function of the theory. As shown in section 5, the two-point function may be written in a Fourier basis. In this basis, Λ is diagonal: ). The right hand side is maximum at ω = 0 (as well as π) and the spectral radius of Λ is κ 2G o (0) 2 . Therefore, for massive fermions and vanishing κ the spectral radius of Λ is κ 2 , whereas for massless fermions we find 1/2d. This way, the leading order approximation is sensible in the vanishing κ and large d limit respectively. In appendix C, we show also that, in the case of strong coupling QCD, where the exact result is known, the effect of higher order terms in Λ in (4.4) vanishes in the large d limit. We turn next to the bosonization of the fermion theory.

Bosonization of fermions
Let A be the matrix with matrix elements: which connects neighbor lattice sites x and y. Then, the fermion action may be split in two terms: with: We integrate fermions by decoupling the quartic term, so that we end up with a quadratic action on fermion fields. If we define ξ(x) matrices with matrix elements: then, the S 1 part of the action is a quadratic form in x, y variables: This action can be written in the diagonal basis of A as: with p being a four vector andξ andÃ representations of ξ and A in the new basis. The exponential function e S 1 is a product of exponential functions. These can be represented using Gaussian integrals: 16) wherec 1 (p),c 2 (p) are integration constants. This way, one finds (modulo an integration constant): 17) whereΣ(p) is a Hermitian matrix. By going back to x representation, the right hand side may be written as an integral with respect to Σ(x, t, t ′ ) fields: Therefore, using (4.12) the action takes the form: Doing the Grassmann integral and using the identity det A = e tr ln A the final expression is: (4.20) Note that A is not invertible and positive definite and this poses a problem in the Gaussian integral representation. However, one may shift A such as to be invertible and positive definite. The cost of doing so is the introduction of another Gaussian field. We have checked that, in this case, the saddle point solution is the same. Therefore, for simplicity we stay with a single Gaussian field and treat A as being formally invertible and positive definite.
In this subsection we have shown that the dual theory is a matrix field theory. Since the action is proportional to N one can solve this theory in the large N approximation. Before doing so we define some useful observables in terms of matrix fields. We start with Green's functions of the theory.

Green's functions
Fermionic Green's function of the theory can be computed by adding the following Grassmann valued source terms in the action (2.19): (4.21) For example, the two point function for t ≥ t ′ , i.e.: may be computed using the double derivative with respect to fermionic sources: Adding fermionic sources in the actionS, equation (4.19), and using the Grassmann integration rules, we get the updated S Σ action: Therefore, the two point function is diagonal in color and d-dimensional space: where: .

(4.26)
Applying the same procedure for the four-point function we find: In the next subsection we turn to the Polyakov loop.

Polyakov loop
A key observable in lattice gauge theory is the Polyakov loop: where the expectation value is evaluated with respect to the effective theory (3.2). In the following we express the expectation value with respect to degrees of freedom of the dual theory. Adding the fermionic source term in the action (2.19): then, the insertion of gauge fields in the path integral may be achieved using the double derivative of the partition function with respect to source fields: where i, j are SU(N) matrix indices and the expectation is evaluated with respect to the action of eq. (2.19). After integration of gauge fields and using the same approximation that led us to the action (4.7) one gets the same action supplemented with source terms: (4.31) The effect of the double derivative in this theory is now the insertion of the following bilinear combination of fermionic fields: where now the expectation value is evaluated in the theory defined by the action of eq. (4.7). Therefore, the Polyakov loop may be written in the form: modulo an irrelevant constant factor and where t o = t N 4 . By means of two point functions this expression may be written in the form: modulo a trivial factor that comes from color summations. In the next subsection, we define one more observable, the fermion-antifermion condensate.

Fermion-antifermion condensate
An important observable in a theory of fermions is the fermion-antifermion condensate: From the discussion in the previous subsection we conclude that: In this study we restrict ourselves in this limited set of observables. In the following, we turn to the solution of the theory in the large N limit.

The large N solution
Since the action, equation (4.20), is proportional to N, one may employ the saddle point solution of the theory. The stationary field should satisfy the necessary first order conditions, which in our case is the system of equations: for any lattice site x of the d-dimensional lattice. The derivative of the first term is taken by expanding the matrix logarithm as a power series on Σ. This way, we obtain the system of equations: These can be inverted to give: Therefore, we get a coupled system of quadratic equations for Green's functions: We are interested in the vanishing κ solution of the theory. In this limit, the second term of the right hand side is small. Given the translation invariance of the time derivative the small κ limit of the solution is also translation invariant, i.e. G(x, t ′ , t) is a function of time separation t ′ − t. Therefore, equations may be written in terms of Fourier transformed Green's functionŝ G(x, ω): for each frequency ω. We solve the system using the solution Ansatz: whereG(x, ω) is a small fluctuation around x-independent solutionG o (ω). Expanding the left hand side of the system (5.5) up to the first order inG(x, ω) and using the fact that (5.5) relates the solution on even and odd sites, we get a scalar quadratic equation: as well as the first order constraint: Substituting the Fourier transformed equation (5.3) in the action (4.20) and writing it in terms ofĜ(x, ω) we get: Expanding the logarithm up to the second order inG −1 oG (x, ω) and substituting the constraint (5.8) we get the effective action in the large N approximation: 10) where translation symmetry ofG(x, ω) on the d dimensional lattice has been used, i.e.
xG (x +μ, ω) 2 = xG (x, ω) 2 . Neglecting cubic order corrections this action describes a free theory of N t bosons with masses given by: The symmetry of the action can be made explicit by writing it in the form of a non-linear sigma model: (5.12) If we defineG(x) matrices with matrix elementsG(x, ω)δ ω,ω ′ , then the global unitary transformations: with U(N t ) denoting the unitary group, leave the spectrum of the theory invariant. In the following subsection we compute the free energy and the mass gap of the theory.

Free energy and mass gap
We compute the solutionG o (ω) of the quadratic equation (5.7). Denoting: we get:G Using equation (5.10) the free energy reads: whereas substitutingG o (ω) in the equation (5.11) the mass spectrum of the theory is given by: The spectrum is bounded from below by the mass gap: In the case of massive fermions and vanishing κ the mass gap M o diverges as m/κ. In the following subsection we study Green's functions of the theory.

Two point function
Using the solution Ansatz (5.6) the leading order approximation of the two point function is: where we have usedγ 5 (x) = ±1 on even/odd sites. SubstitutingG o (ω) from equation (5.15) and expanding the square root in κ for massive fermions we get: In the following we concentrate on even sites two-point function, i.e.γ 5 (x) = 1 and compute the first term of the right hand side. The odd sites expression is then the negative even sites result with the formal substitution m → −m. The time domain two-point function may be computed using the inverse discrete Fourier transform: plus O(κ 2 ) terms. The leading term X + (t) is thus the propagator of free fermions in 0+1 dimensions. The latter may be obtained by solving the linear system of equations: with antiperiodic boundary conditions. Both ways we find: 4 The right hand side of (5.21) may be computed in two steps: first, one computes the N t → ∞ expression, i.e. X is the Heaviside function. Then, the finite N t result is obtanied by evaluating the infinite sum X + (t) = +∞ m=−∞ (−1) |m| X (∞) with sinh E = m, C = cosh E cosh(EN t /2) and N t even. In case N t is odd the above result should be replaced by the expression: The different behavior of the two point function at even and odd times is a manifestation of fermion doubling on the lattice, i.e. the presence of two poles of 1/(m + i sin ω) at iE and π − iE corresponding to the same energy E.

Polyakov loop
We begin by expressing the time domain trace of (4.34) in the frequency domain: In the large N approximation the right hand side factorizes and we find: since for N 4 even phase factors cancel due to equal number of even and odd sites in the product. If N 4 is odd then only one phase factor remains. Due to the π-periodicity ofG o (ω) the Polyakov loop may be written as a sum over positive frequency terms: 5 We evaluate the sum in the large N 4 limit, in which case the largest term dominates. Note thatG o (ω) is maximum at ω = π/N t . In the large N t limit π/N t is close to zero and we evaluate: This way, the free energy F o of the static charge is: where we have restored the lattice spacing a. At leading order in κ and m = 1 we havẽ G o (0) = 1 − 2dκ 2 + O(κ 4 ) and therefore: which vanishes in the limit κ → 0. The result does not change if we useG(π/N t ) instead of G(0) as well as if we add the next largest term in the sum over frequencies in (5.27).
Using F o we can compute the renormalization group the beta function: It is negative and vanishes linearly with the coupling constant. The theory is thus asymptotically free. However, its ultraviolet behavior is different from the standard Yang-Mills theory. The same conclusion may be drawn if we had computed Wilson loops. In this case we would find a perimeter law.

Fermion-antifermion condensate
Using its definition, equation (4.36), and substituting the two point function in the leading order approximation of the large N solution, the condensate is: where we have used again the π-periodicity ofG o (ω). Expanding the right hand side in κ and taking m = 1 we have: For large N t the first term is a lattice sum equal to N t / √ 2 6 and therefore: This way, in the continuum limit, the theory is characterized by a non-zero value of the fermion-antifermion condensate.

Synthesis
The solution of the theory in the large N approximation shares distinctive properties with the standard Yang-Mills theory like asymptotic freedom and color confinement. Note however, that the beta function of our effective Yang-Mills theory is different from that of the standard Yang-Mills theory. In the next section we probe the theory in the massless limit.

Asymptotic safe QCD
We have seen that the effective Yang-Mills theory is local and asymptotically free. In this section we compute the solution in the case of massless fermions, or almost massless fermions with mass m. In subsection 3.2 we learned that the effective theory of our model is QCD with a large number of flavors. The large N solution is found following the same steps as in the previous section. The mass spectrum of the theory is given by equation (5.17) with: However, the lightest mass: vanishes for exactly massless fermions. Therefore, the theory is gapless. Note however, that the number of light fermions is large, as can be seen from equation (3.12). If there are N f such modes, then the effective action of the theory, equation (5.12), may be split in two pieces corresponding to light and heavy modes: 3) with Ω l and Ω h being the set of light and heavy mode frequencies. If we define light mode matricesG l (x) with matrix elements: and take light modes to be massless, the action is symmetric with respect to global U(N f ) L × U(N f ) R chiral transformations of light modes: As shown below, the fermion-antifermion condensate of the full theory is nonzero due to light modes alone. Therefore, the chiral symmetry of the light modes is spontaneously broken to U(N f ). Taking the limit of vanishing m in equation (6.2) the mass of the N 2 f Goldstone modes is: a result which is expected from the chiral perturbation theory [22]. Therefore, the solution shows that the low lying spectrum behaves as in QCD. Note however that N f is large in our case.

Two point function
We compute the two point function using again the solution Ansatz (5.6) in the leading order approximation. The massless two point function is then: where the plus/minus subscript corresponds to even/odd sites of the d dimensional lattice. SubstitutingG o (ω) from equation (5.15) we get: In continuous time, this is the same as the two point function of the q = 2 SYK model, i.e.: where J is the coupling constant of the SYK model, with the identification J = √ 2dκ. The time domain two point function is given by Maldacena and Stanford [15]: 2dκ|t| sin θ , (6.10) which for large time separations is: This is a power law decay as opposed to exponential decay in the case of massive fermions.

Asymptotic safety
Following the same steps as in the case of massive fermions we compute the Polyakov loop, equation (5.27), in the large N 4 limit. The free energy of the static charge is thus evaluated at zero frequency using the same formula (5.29). The result in the massless case is: Using F o , the beta function of the theory is: It is zero at κ = 0 and κ c = 1/ √ 2d, positive for κ ∈ (0, κ c ) and negative for κ > κ c . Hence, the theory has an ultraviolet fixed point at κ c . Solving the renormalization group equation: withm being an integration constant, the correlation length of the theory is defined by: At the critical point, it diverges according to the law: i.e. the theory shares the same critical exponent with two-dimensional Ising model. The theory has a continuum limit at a critical value of the coupling constant. Note that in the limit d → ∞ the theory becomes asymptotically free.
There are some consequences of the critical theory. For exmaple, at critical κ the mass of the Goldstone boson, equation (6.6) becomes: which again vanishes as expected from chiral perturbation theory. Another consequence is that the length of the extra dimension is finite in the continuum limit. However, due to Therefore, both theories considered in this paper, for massive and massless fermions, satisfy the large N t assumption which has been used by us from the beginning.

Fermion-antifermion condensate
Splitting the sum in the equation (5.32) in to light and heavy frequency contributions we find: (6.14) Therefore, in the massless limit we get: This way, in the continuum limit, the theory has a non-zero chiral condensate. At critical κ its value is independent of the coupling constant and the number of dimensions.
1) where G(t, t ′ ) and Σ(t, t ′ ) are bilocal fields defined on a time lattice, N is a large positive integer and J is the coupling constant of the theory. In this subsection we relate our model, i.e. equation (4.20): 2) to the q = 2 SYK model. The relation is established in the massless case under the following conditions: i) the large N limit; ii) the coupling constant relationship J 2 = 2dκ 2 .
The saddle point solutionG o (ω), which is x independent, suggests that we may approximate the matrix A by a constant matrix, i.e. A = 2d. This way, the action (4.20) decouples completely in x-space: This is an ideal gas of pairs of one-matrix theories: corresponding to even/odd sites of the d dimensional lattice. Picking even sites only and setting m = 0 we end up with model: (7.5) One-matrix models have been studied in the past as a non-perturbative formulation of the two-dimensional gravity, see for example [23]. In general, such models do not lead to black hole formation [24]. As it was shown in the previous section, the theory with massless 7 See for example equation A.9 of Gross and Rosenhaus paper [14].
-25 -fermions shares the same two-point function with the action of q = 2 SYK model. This is not a coincidence. The action (7.5) may be written in terms of an additional Hermitian matrix G(t, t ′ ): (7.6) which can be shown by integrating e S Σ,G,1 with respect to individual matrix elements G(t, t ′ ) using Gaussian integrals (4.15) and (4.16). Rescaling the matrix G → −iG as well as the matrix Σ → −Σ we get twice the action of the q = 2 SYK model on a time lattice: 7) where 2dκ 2 = J 2 is identified with the square the coupling constant of the SYK model. Therefore, the large N asymptotic safe QCD may be described as an ideal gas of q = 2 SYK models on each site of a d dimensional lattice. Note that since d and κ are related in continuum limit such that κ c = 1/ √ 2d we have J 2 c = 1. Therefore, the relationship (3.4), i.e. N t = 4κ −5 , has no influence on the magnitude of J. It merely tells that at critical κ the length of the extra dimension N (c) t = 4(2d) 5/2 is large. Therefore, the continuum limit of the asymptotic safe QCD corresponds to the low temperature limit of the q = 2 SYK model with J 2 = 1.
Note that for q = 2 the SYK model is not chaotic, whereas for q = 4 it saturates the chaos bound. Since we used the leading order approximation of the F function (see equation (4.6)), the effect of order q = 4 terms or higher is expected to be present in the effective action. The extent to which these terms alter the chaotic behavior of our theory remains unclear without a proper calculation. However, the q = 2 SYK model originates from a quadratic Hamiltonian with disorder couplings. Magan has shown that a generic model of quadratic fermions with random couplings satisfies the Eigenstate Thermalization Hypothesis [25].

Summary and discussion
In this paper, we have formulated and studied a lattice theory of fermions beyond the Standard Model. Integration of fermions yields a lattice gauge theory which is expressed in terms of Wilson loops of growing sizes. This premise is interesting alone since the effective theory is a local Yang-Mills theory in the limit of vanishing coupling constant. On the other hand the fermion theory may be integrated with the help of known one-link integrals. The remaining effective theory of fermions is then bosonized with the help of Hermitian matrices. The dual theory obtained this way is solved in the large N limit. The solvability of the theory is a distinctive property of the model. The solution shows that the model shares qualitative properties of strong interactions like asymptotic freedom, color confinement and a spectral gap. However, the renormalization group beta function of the theory vanishes linearly with the coupling constant as opposed to the cubic law of the standard Yang-Mills theory [1,2]. Nonetheless, the main result of the paper is that a local Yang-Mills theory exists which is non-perturbatively solvable in contrast to the present status of an unknown similar solution to the standard Yang-Mills theory in four dimensions.
We have studied also the model with massless fermions. In this case the effective theory of Yang-Mills fields is QCD with a large number of light flavors. Its solution is ultraviolet complete at a non-zero critical coupling constant where the theory is scale invariant. The light modes of the theory are shown to be chirally symmetric, a symmetry which is spontaneously broken. The theory shares the same two point function with the q = 2 SYK model in the leading order of large N approximation and the leading order approximation of our bosonization approach.
In the appendix A we have shown that one may use complex N × N matrices as disorder couplings instead of SU(N) couplings. Within the approximation made for the bosonization of fermions in case of the SU(N) disorder and the scaling relation assumed between N and N t (see appendix A) we show that SU(N) and Gaussian disorder give closely related effective theories.
The results of this paper show that the model has a rich structure, which we intend to study further in the future. The next step is to treat the bosonization of fermions exactly. In order to probe the theory further we would like to compute more physical quantities such as the low lying meson spectrum. The computation of fermion-antifermion potential would reveal the nature of the interactions of the theory. We would like to study further the connection to the SYK model with the intention to find whether the theory has a gravity dual. Finally, we would like also to simulate the theory on the computer.

A Gaussian disorder
The SU(N) disordered fields are a special case of more general model of lattice fermions coupled to general N × N complex matrices Φ µ (x) attached at each directed link (x, x +μ) on the lattice. Adopting the same notations as in subsection (2.1) the Hamiltonian operator of the theory is: x,a,b,µγ where matrix elements of disorder fields, i.e. Φ µ (x) ab , a, b = 1, 2, . . . , N are distributed according to density: with C being a normalization constant. The theory may be formally described by the same Hamiltonian kernel h as in subsection (2.1) with SU(N) fields substituted by Gaussian fields.
In terms of the Grassmann valued fermion field ψ(x, t), with t labeling points in the extra dimension, the model is defined by the action: whereas the partition function of the theory is: Gaussian fields can be integrated using the formula: where α is a positive real number and β, γ are complex numbers. The effective action of the theory which remains after integration of Gaussian fields is: This is precisely equation (4.7), which is used as an approximation of the full fermion theory (4.3) in the case of SU(N) disorder. This shows that SU(N) and Gaussian disorder theories are closely related. In the following subsection we elaborate further this relationship.

A.1 Embedded gauge fields
In this subsection we identify embedded gauge fields within Gaussian fields and show that the large N t theory is effectively a local Yang-Mills theory. We begin by giving the expression of the effective action which remains after fermions are integrated out. Following the same steps as in subsection (3.1) and taking into the account the Gaussian measure (A.2) the effective action of the theory is: Note that the trace in the right hand side of (A.5) is taken in the tensor product space of the lattice sites and N × N matrices. Using the polar decomposition of Gaussian fields: gauge fields are identified by the U(N) factor U µ (x), where φ µ (x) are positive definite Hermitian matrices. In the following we keep the gauge fields fixed and find the saddle point action in the large N t limit using the solution Ansatz: where ϕ is a real value. A more general Ansatz would include terms which give O(1/N t ) contributions to the saddle point action. Since we are interested in the leading contributions to the effective action we stay with the above solution Ansatz. Note also that both N and N t are large and we relate them by the κ dependent factor: Using this definition and the solution Ansatz the action is a function of a single variable ϕ: where h o is the fermion matrix in the background of the gauge field U µ (x). The saddle point equation S ′ (ϕ) = 0 yields the nontrivial solution ϕ o given implicitly by the equation: Note that the solution ϕ o is gauge field dependent. As shown in the next subsection, an explicit solution may be computed in the limit of vanishing κ. Alternatively, one may assign to the action (A.10) an approximate as well as gauge field independent solution. Such a solution corresponds to an approximate saddle point. This way, one may proceed as in subsection (3.1) and obtain a local theory of Yang-Mills fields. Therefore, Gaussian and U(N) disorder theories are related at the approximate saddle point of the Gaussian disordered action in the large N t limit. In the next subsection we give an example of a precise relationship.

A.2 Saddle point Yang-Mills theory
In this subsection we compute the saddle point action of the theory with Gaussian disorder fields in the limit of vanishing κ. Our starting point is the saddle point equation (A.11).
Making the following Ansatz for the left hand side: with a 2 , a 4 being real constants, and matching it to the right hand side expansion of (A.11) in κ we find: Note that we have neglected the higher powers of the expansion of the right hand side since we seek the limit of vanishing κ. We have also the freedom to select a small value of a 4 in order to control the matching error of the expansion. Substitution of a 2 and ϕ 2 o to the action (A.10) give the effective theory of Yang-Mills fields: Since the leading term does not look like the standard plaquette action of Wilson we expand Tr h 4 o in terms of gauge fields and get: For smooth gauge fields close to continuum limit, the plaquette terms TrU µ U ν U * µ U * ν are close to one. This property allows us to write the right hand side as a geometric series of Nd with the leading term: If we insist to maintain the scaling of N t to κ as defined in (3.4) then the leading terms of this theory and the one of (3.5) are the same provided we select a 4 = 3/4. This calculation shows that close to continuum limit the saddle point action of the theory with Gaussian fields yields a similar Yang-Mills leading term as in the case of the theory with SU(N) disorder fields. Note that the form of the r(κ) Ansatz (A.12) is crucial to arrive to this conclusion. If there is no relation between N and N t the theory may not be local.

B Group integration
This section is written to make the paper self contained. We begin first with some integration rules in the unitary groups.

B.1 Group integration rules
Unitary group integration rules used in this paper rely on the Haar measure. These rules are known and we point the reader to the paper of Creutz, reference [16], 8 for a detailed account.
Here we would like to give a few useful results. For example, invariant group integration arguments lead to the conclusion: One can extend this result for the product of n matrix elements as long as n = N. In the case n = N the integral is nonvanishing if indices lead to a group invariant quantity such as the determinant of U, 9 i.e.: where a 1 , a 2 , . . . , a N is a permutation of 1, 2, . . . , N and ǫ a 1 a 2 ...a N is the rank N totally antisymmetric tensor. Indeed, if we multiply both sides by ǫ a 1 a 2 ...a N and sum over all permumations a 1 , a 2 , . . . , a N we get an identity. Another useful integral is: One can be convinced about the normalization and the δ ad factor by taking b = c and summing both sides over b. The same argument justifies the factor δ bc . This integral is also derived at the end of the next subsection applying the results obtained therein.

B.2 Computation of one-link integrals
In this subsection we deal with the computation of one-link integrals of the type: where a, b indices run from 1 to N, the t index runs from 1 to N t , the trace is taken in t-space and dU is the Haar measure of the U(N) group. Using invariance properties of the Haar measure, the integral depends only on gauge invariant quantities: On the other hand, the integral depends also on bilinear Grassmann sums: since they are invariant with respect to invertible N t × N t matrix transformations. Therefore, if we define the N t × N t matrix: the integral is a function of t-space traces of this matrix: tr Λ , tr Λ 2 , . . . , tr Λ Nt , (B.8) i.e. it can be written in the form: where F is a matrix valued function defined by its power series expansion. Taking the derivative of both sides with respect toψ t 1 b and then ψ t b in this order, multiplying byψ t 2 a ψ t a , summing over a, b and t we find: The left hand side may be written in terms of derivatives of tr F (Λ) with respect to matrix elements of Λ. This way, we obtain a system of N 2 t coupled second order differential equations for trF (Λ): supplemented by suitable boundary conditions. For example, the set: guarantees that the solution is regular around zero. Using the definition of Λ, equation (B.7) the system to be solved is: This system may be written in matrix notations in the form: Since we are interested in the large N limit we drop the third term and find: This is an algebraic quadratic matrix equation for Λ ∂ trF ∂Λ . Its solution poses no problem and is given by: where the second condition in equations (B.12) is taken into account and the matrix square root is defined in terms of its power series expansion. A solution for the matrix F (Λ) that satisfies this equation as well as the condition F (0) = 0 is: In this case too, the matrix logarithm is defined in terms of its power series expansion. 10 As an application let us derive the group integral (B.3). Taking the fourth derivative of both sides of the integral (B.9) with respect toψ, ψ,χ, χ Grassmann variables in this order one finds: Substituting F (Λ) in the left hand side using the result, equation (B.17), one finds: Comparing right hand sides of last two equations, equation (B.3) is thus established. Other interesting integrals can be computed using this method.

B.3 Relation to other work
Our derivation can be related to the one of reference [20]. In this reference the derivative of tr F is taken with respect to invariant traces: (B.20) The resulting system of differential equations is related to ours, equation (B.14), if the derivative of tr F with respect to matrix elements of Λ is defined by the expression: In order to find the solution, reference [20] diagonalizes the matrix Λ, whereas reference [21] relies on the "strong coupling" solution of Brezin and Gross [18], their solution being found also by diagonalizing Λ.
C QCD at strong coupling The large N limit of QCD at strong coupling has been studied in the past. Notable references are Kluberg-Stern et. al. [20] and Kawamoto and Smit [21]. The aim of this appendix is to derive the main results of QCD at strong coupling using the bosonization approach employed in this paper. The action is given by equation (2.19) but without time derivative term, i.e.: where κ is now fixed at the value of 1/2. Everything else being the same, integrating gauge fields as in subsection 4.1 we get: where F (.) is defined by the expression (see (B.17)): In this paper we approximate the right hand side with the leading order result F (Λ) = −Λ + O(Λ 2 ). It is this approximation that will be tested in the case of strong coupling QCD. Following the same steps as in section 4 the final expression of the bosonic effective action is (see equation (4.20)): For future references on QCD at strong coupling we make this section self contained. Therefore, a few steps of the large N solution of section 5 will be repeated here. In order to find the large N solution we find first the field that makes the action stationary and then compute an effective action by computing fluctuations around such a solution. The stationary field should satisfy the necessary first order conditions, which in our case is the system of equations: where S Σ is the bosonic action, equation (C.4). The derivative of the first term is taken by expanding the matrix logarithm as a power series on Σ. This way, we obtain the system of equations: Denoting by G(x) = 1/[m f + Σ(x)] the N t × N t matrix, equations take the form: This is a system of matrix valued quadratic equations which is solved using the Ansatz: whereG(x, t, t ′ ) is a small fluctuation field around the uniform solution G o δ t,t ′ . Expanding the left hand side of (C.7) up to the first order in power series of the matrix G −1 oG (x), we get a quadratic equation for the matrix G o : with the solution: as well as the first order constraint in the fluctuating fieldG(x, t, t ′ ): Substituting the Ansatz, equation (C.8), into the effective action (C.4), using equations (C.6), (C.7) and the first order constraint (C.11), as well as expanding the logarithm up to the second order in powers of the G −1 o G(x) matrix we get: where the translation invariance on the lattice, i.e. the identity xG (x+μ, t, t ′ ) 2 = xG (x, t, t ′ ) 2 has also been used. From this expression one can infer the free energy of the theory: On the other hand, equation (C.12) can be written in the form: with the trace taken in t-space and where we have denoted: Neglecting cubic order corrections, the resulting effective theory (C.14) is a non-linear sigma model of massive bosons with mass M given by equation (C.15). For vanishing m f , the mass squared vanishes linearly with m f : For vanishing quark mass the theory has an exact global U(N t ) L × U(N t ) R symmetry, i.e. the effective action is invariant with respect to global transformations: The chiral condensate of the theory does not vanish and therefore the chiral symmetry (C.17) is spontaneously broken to U(N t ). In order to see this, we compute the chiral condensate, which is defined by expressions: (C.18) Using the result (C.13) for the free energy and substituting the solution G o , equation (C.10), we get: Therefore, the chiral symmetry of QCD at strong coupling is broken spontaneously to U(N t ) group, whereas the spectrum has N 2 t Goldstone bosons. We note that our results rely on the leading order approximation F (Λ) = −Λ + O(Λ 2 ). Our condensate and boson mass differ with the exact results of references [20,21] by a O(1) prefactor at vanishing fermion mass. This small difference does not alter the physical picture and thus justifies our approximation. In order to make these differences transparent, we summarize below the approach and results obtained in these references.

C.2 Relation to other work
If we want our results to be exact, we may follow the bosonization approach of references [20,21]. Another possibility is to use the Lagrange multiplier method as shown in a series of papers related to the SYK model, see for example Gross and Rosenhaus [14], as well as Kitaev and Suh [13].
Kluberg-Stern et. al. paper starting point is the observation that integration of functions which depend on Grassmann pairsη t η t ′ , t, t ′ = 1, . . . , k can be cast in the form of the following matrix determinant:  . (C.20) The proof can be established by noting that only the order k term of Grassmann pairs power series expansion of f contributes in the integral. If Grassmann pairs of the right hand side are formally substituted by real valued matrix elements σ tt ′ = 0, t, t ′ = 1, . . . , k we get the identity: Substituting the Fourier representation: in the right hand side we get: The Fourier transformed functionf can be written in terms of the original function using the inverse Fourier transform and the final expression is: (C.24) If f depends on Grassmann sums N a=1η t a η t ′ a , t, t ′ = 1, . . . , k, as it is the case in our application, the Berezin integral, equation (C.21), gives: This way, the fermion theory can be cast in the form of the bosonic theory with an action S which depends on matrix elements λ tt ′ , σ tt ′ , t, t ′ = 1, 2, . . . , k: − S(λ, σ) = tr (N ln λ + iλ T σ) + ln f (σ) . (C.27) Using this technique and adopting the normalization of fermion bilinears as in reference [20], i.e.ψψ → iψψ, one can show that the large N solution of the fermion theory (C.2) is given at the saddle point σ + δ tt ′ with: whereas the chiral condensate is: Now let us turn to the approach followed by Kawamoto and Smit [21]. These authors use the bosonization of the fermion determinant based on the identity: where dU is the Haar measure on the U(k) group, R is a Hermitian and positive definite k ×k matrix and c k N are defined by the expression: c k N = 0! 1! · · · (k − 1)! N! (N + 1)! · · · (N + k − 1)! .   The proof of the identity (C.30) is given in the appendix of the reference [21]. However, the basic idea can be understood in the case N = 1. Indeed, setting R = 1 and factoring U = e iϕ V , where V is a SU(k) matrix one has: dU e trJU det U = 1 2π Expanding the second exponential and integrating over ϕ only the power of order k gives a nonvanishing result. Therefore, one gets: where we have used the group integration formula: whereas the free energy of the theory is: This way, the chiral condensate is given by the expression: Therefore, whichever bosonization method is used one gets identical results.
Finally, let us compare these results with ours in the limit limit m f → 0. Our saddle point solution, equation (C.10), evaluated at κ = 1/2, as well as the condensate, equation (C.19), read: (C.39) Therefore, our approximation misses the prefactor 1 − 1 2d , which is of the order O(1) even at d = 2. For large d, the exact treatment of the function F , equation (C.3), yields the same result as ours.