Solvable Quantum Grassmann Matrices

We explore systems with a large number of fermionic degrees of freedom subject to non-local interactions. We study both vector and matrix-like models with quartic interactions. The exact thermal partition function is expressed in terms of an ordinary bosonic integral, which has an eigenvalue repulsion term in the matrix case. We calculate real time correlations at finite temperature and analyze the thermal phase structure. When possible, calculations are performed in both the original Hilbert space as well as the bosonic picture, and the exact map between the two is explained. At large $N$, there is a phase transition to a highly entropic high temperature phase from a low temperature low entropy phase. Thermal two-point functions decay in time in the high temperature phase.


Introduction
In this paper we are interested in the physics of a large number of non-locally interacting fermionic degrees of freedom. We will study quantum mechanical fermions with a vector-like index structure as well as a matrix-like index structure. Part of the motivation comes from recent investigation of similar systems [1,2,3,4,5,6,7,8] which may play a useful role in understanding certain problems in black hole physics.
Another motivation comes from recent investigations on the emergence of bosonic matrix models from discrete systems [9,10,11,12,13]. Moreover, we are interested in how the original fermionic degrees of freedom and Hilbert space might be encoded (if at all) in those of the bosonic matrix. 1 Our goal is to study tractable systems with rich interactions. Unlike the purely fermionic models of [2], our models do not have any quenched disorder (see also [14,15,16,17]). Instead, their complexity stems from the matrix-like interactions. Given a matrix structure, we might hope, eventually, to relate our systems to some gravitational description for a new class of models.
Here, we take some small steps in this direction by carefully analyzing and solving such systems with quartic interactions. Our approach is to express the systems in terms of auxiliary bosonic variables allowing us to use standard large N techniques.
When possible, we carry forward the analysis in both the fermionic Hilbert space picture as well as the bosonic path integral picture and map physical questions from one picture to another. The bosonic systems exhibit a (0 + 1)-dimensional emergent gauge field and the corresponding (broken) gauge symmetry is crucial in our ability to solve the systems. Part of our treatment is in some regard a (0 + 1)-dimensional analogue of the analysis in [18,19]. Fermionic correlation functions are related to the calculation of certain Wilson line operators of the gauge field. At the end, for both matrix and vector models, we are able to express the exact finite N thermal partition function as an ordinary integral. In the matrix case, this integral is a matrix-like integral with a modified Vandermonde term that consists of the sinh or sin of the difference in eigenvalues rather than the difference eigenvalues themselves.
Such eigenvalue integrals appear in the study of unitary random matrix models as well as Chern-Simons theories. At large N the systems exhibit a non-trivial phase structure which is naturally characterized, in the matrix case, by the connectivity of the eigenvalue distribution, somewhat similar to the situation encountered in [12].
The structure of the paper goes as follows: we begin by analyzing the vector model in the early sections. We show how the fermionic thermal partition function, initially expressed as a path integral over Grassmann valued functions of Euclidean time, can be expressed in terms of an ordinary bosonic integral over a single real variable. We analyze both the thermal phase structure and fermionic correlation functions. The latter are expressed in terms of the expectation values of non-local in time variables, resembling Wilson line operators of an emergent gauge field. We show, at large N , a transition from a low entropy phase to one with entropy extensive in N . These results are generalized to the matrix case in the latter sections, where the structure is significantly more intricate. The thermal partition function now reduces to a matrix integral, which as mentioned, contains a modified Vandermonde interaction among the eigenvalues. We end discussing the thermal phase structure and thermal correlation functions of the matrix model at large N . At large N , we find that the real time two-point function decays significantly faster in the matrix case than the vector case. In appendices A and B we discuss certain generalizations of the models studied in the main body.

Fermionic vector model
In this section we consider N complex Grassmann degrees of freedom {ψ I ,ψ I } interacting via a quartic interaction. The index I = 1, 2, . . . , N is a U (N ) index, under which the ψ I transform in the fundamental and theψ I in the anti-fundamental representation.
The thermal partition function of our model is given by: We have a periodic Euclidean time coordinate τ ∼ τ + β, such that our integration variables obey ψ I (τ + β) = −ψ I (τ ). The fermion variables are dimensionless and the parameter γ is a positive number with units of τ . (In appendix B we carry over our results for the vector model to the case with negative γ.) Using a Hubbard-Stratanovich transformation the Euclidean action can be written as with λ(τ ) a periodic function of τ with units τ −1 . The only dimensionless quantity, other than N , is γ/β. It grows with increasing temperature. We set β = 1 unless otherwise stated.

Exact evaluation of path integral
In this subsection we evaluate the fermionic path integral exactly. Going back to (2.2), evaluation of the fermionic path integral yields: The functional determinant has a local invariance [11]: This can be viewed as a time-reparameterization of an einbein λ(τ ), and it is an exact symmetry of a system of non-interacting fermions. Due to the above invariance, the functional determinant will only depend on the Matsubara zero-mode of λ(τ ), namely where the Green function G n,m contains the diagonal component of the λ n,m : andλ has Toeplitz form:λ We can expand the logarithm on the right hand side of (2.5) in a matrix Taylor expansion: (2.8) The zeroth order piece, using the standard Euler infinite product formula, can be written as 2 The first non-trivial contribution inλ to (2.8) is quadratic and involves n∈Z 1 (2πi(n + 1/2) + λ 0 ) 10) The vanishing of the above expression implies that no derivatives of λ(τ ) are generated to leading order. It is not hard to check that this behavior continues to higher orders.
Thus, the path integral (2.3) in a thermal frequency basis becomes Performing the integration over the non-zero modes, we arrive at: where we have fixed the normalization constant in (2.11) by demanding that in the infinite temperature limit we get the Hilbert space dimension Z = 2 N . We see that the partition function can be reduced to a simple integral over λ 0 (the unique quantity invariant under (2.4) that can be constructed out of λ(τ )). The integral over λ 0 can be explicitly done. Reinstating β, we find: 13) where C N n are the binomial coefficients. From Z[β] we can read off the spectrum and degeneracies: (2.14) Note that the spectrum is symmetric under n → (N − n). 2 The infinite constant is absorbed in the normalization factor N in (2.3).

Hilbert space picture
Can we understand the above result from the point of view of the Hilbert space? Upon quantizing the system, we impose anti-commutation relations {ψ I , ψ J } = δ IJ .
If we define an empty state |0 as annihilated by all the operators ψ I , then the full set of states is given by acting with any amount ofψ I on |0 . This gives 2 N states.
The Hamiltonian of the system is given by: where we have fixed the normal ordering constants such that the spectrum matches that of (2.14). The spectrum is non-positive for γ > 0. The doubly degenerate ground state has energy E g = −N/(16γ) . Since the operatorN = ψ I ψ I commutes with the Hamiltonian, we can organize the spectrum in terms of eigenstates of thê N operator, which counts the number ofψ I hitting |0 . For exampleNψ 1ψ2 |0 = 2ψ 1ψ2 |0 . Thus, states with nψ's acting on |0 have an energy E n given in (2.14) and a degeneracy given by the binomial coefficients d n = C N n . As expected, if we sum over all the degeneracies with find: n C N n = 2 N . At large N the spectrum peaks sharply about n = N/2, with d N/2 ∼ 2 N +1 / √ 2πN , consisting of states with E N/2 = 0.

Thermodynamic Properties
From the thermal partition function we can study various thermodynamic properties.
At very low temperatures the entropy is O(1) while the energy goes as O(N ), and the free energy F = −T log Z = E −T S is dominated by the energy piece. As we increase the temperature, the energy and entropy contributions begin to compete since the number of states begins to grow exponentially in N . More precisely, at large N the degeneracies behave as d n ≈ C N N/2 e −(N −2n) 2 /(2N ) . When d n = e Sn ≈ e βEn , which at large N gives β ≈ 8γ, we expect a transition. Above this temperature, the entropy become O(N ) and dominates over the energy. The transition becomes increasingly sharp as we increase N . We show a numerical example in figure 1. Note that in the high temperature phase, the entropy of the system is extensive in the number of degrees of freedom.
as a function of β. We have taken N = 500 and γ = 1. Notice the transition occurring near β = 8γ.

Correlation functions of the vector model
In this section we discuss the correlation functions of the model in both real and Euclidean time. We do this both in the path integral picture and the corresponding Hilbert space picture.

Fermionic Hilbert space picture
To compute the thermal correlator in the Fermionic Fock space picture, recall that we can organize the Hilbert space in terms of number eigenstates: |n, I n =ψ I 1 . . .ψ In |0 , m, I|n, I = δ m,n δ I,I . The real time two-point function in the thermal ensemble is explicitly given by: The (N − 1) in C N −1 n comes from the reduction in the number of states in |n, I n when hit with the operatorψ B . We note that the Green function (3.2) satisfies At low temperatures, the two-point function oscillates with a frequency given by we expect to see recurrent patterns in the two-point function for time separations of the order t ∼ N γ. In figure 2 we display a plot of G β (t) exhibiting this behavior.

Path integral expression
The generating function for Euclidean fermion correlation functions is given by: Integrating out the fermions we get Consequently, the exact Euclidean time two-point function is given by: The differential operator: obeys the following equation: For 0 < τ, τ < β, we find: At this point, we must fix the remaining constant. This follows from the fermionic nature of the correlator under a shift in β, i.e. G(τ 1 + β, τ 2 ) = −G(τ 1 , τ 2 ). We thus have: Notice that G(τ, τ ) is now only invariant under the transformations (2.4) that leave the end points unchanged. Recalling the interpretation of λ(τ ) as an einbein, we can view G(τ, τ ) as a gravitational Wilson line.
Since G(τ, τ ) can be expressed as the exponential of a local integral of λ(τ ), we are now in a position to evaluate the path integral. Upon evaluation of the functional determinant, a cancellation occurs between c + and one of the powers of the cosh in (2.12), and we can express the path integral as: Performing the Gaussian integrals pointwise, and reinstating β we find for 0 <τ < β: The normalization constant N has been fixed by imposing G AB β (0) = 1/2. This agrees precisely with (3.2) upon Wick rotating to Euclidean time t → −iτ .

Large N approximation
We would like to see how much of the exact structure previously uncovered is contained in a large N approximation. It is convenient to express the partition function in terms of thermal Fourier modes and reinstate β. From (2.9) and (2.11) we get: (3.13) The large N saddle point equations for λ(τ ) are: amounting to a constant solution for λ(τ ). For low temperature one finds three saddles, which to leading order in the large β limit are: λ 0 = 0 and λ 0 = ±1/4γ, the one at the origin being subdominant. As we increase the temperature and reach β = 8γ the λ 0 = 0 saddles coalesce to the origin. Recall that β = 8γ was the temperature at which the free energy exhibited a thermal transition. For large temperatures, β < 8γ, a single saddle is found corresponding to λ(τ ) = 0.

Low temperatures
In the Hilbert space picture, the dominant behavior in the β → ∞ limit comes from the vacuum and its first excited state. From (3.2) we get in the large N limit which shows a single oscillatory behavior in time with frequency given by the energy difference of the first excited state with respect to the vacuum. In the large N limit this is given by ∆E = 1/4γ .
This behavior is recovered, in the path integral description, from the saddles at with ω p = 2π(p + 1/2)/β and p, q ∈ Z. We can expand the Green function as: To leading order in N we therefore keep the 1 in (3.18). In the small temperature limit we transform back to τ replacing p → β dω and we also have log cosh βλ 0 2 ≈ β|λ 0 | 2 . Evaluating (3.17) at the saddles we get which coincides with (3.16) when Wick rotating to real time τ → it. Notice that τ > 0 (τ < 0) implies that only the λ 0 = +1/4γ (λ 0 = −1/4γ) saddle contributes to the ω contour integral.

High temperatures
Consider first the Hilbert space picture. The high temperature limit of (3.2) is dominated by the n ≈ N/2 terms due to the high degeneracy of states. At large N we can and defining the variable x = −(N − (2n + 1)))/N we get: Again, the normalization constant c was fixed by demanding that G AB β (0) = 1/2. In the second line, we have kept only the leading in N terms in the exponent. Interestingly, the correlation functions decay to exponentially small values after a time of order t ∼ N γ(γ − β/8). This is in spite of the correlator not exhibiting an initial exponential decay of the form e −t/β , characteristic of thermal systems. As mentioned below (3.3), we expect the approximation to fail and find recurrences when all the oscillating factors near n ≈ N/2 are in phase. This happens for t ∼ 4πN γ. So, the recurrences in the vector model occur much more frequently than for a strongly coupled (chaotic) system [27], where they might be separated by exponential in N (or even super-exponential) time scales. These features suggest that the system has a flavor of integrability. The approximate result (3.20) agrees well with the exact answer (3.2) for times t 4πN γ.
We now consider the path integral picture. From (3.6) we have The auxiliary parameter, a, is introduced to make the exponent local in τ . The approximation in the second line corresponds to taking the β → 0 limit. In this limit λ 0 1 and we can approximate the log cosh λ 0 /2 by a quadratic approximation. We can now perform the Gaussian integral (3.22) as was done in (3.12), getting where again we have fixed the normalization such that G AB β (0) = 1/2 at τ = 0. This coincides with (3.20) upon Wick rotating. Note that in the large N approximation used above, we have gone being the leading saddle point approximation for which λ(τ ) = 0, leading to a constant correlation function in τ . The τ -dependence in (3.22) arises from the next to leading (Gaussian) correction about the large N saddle point.
As a final note, we should mention that at β = 8γ, the correlator exhibits a transition between the high temperature decaying correlator and the low temperature oscillatory one. At large N and β = 8γ, we must consider the log cosh λ 0 /2 term beyond the quadratic approximation. For instance, keeping only the quartic term we are able to approximate the correlator at β = 8γ. In real time, it takes the form: (3.23) The above agrees well with numerics at large N , as seen in figure 3. It falls to small values after t ∼ N 1/4 β, which is parametrically faster than the high temperature case.
Our large N analysis is generically not sensitive to the presence of recurrences. Indeed, for large times (3.22) will receive significant corrections from the terms in (3.18) that were discarded. This requires us to go beyond the saddle point approximation.
These are responsible in restoring the fine structure of the correlation function. 3

Alternative view of the low energy sector
Here, we provide an alternative view to the fermion determinant independence of the non-zero modes of λ(τ ) (2.4) and describe the theory's effective dynamics at low temperature. This approach will be useful when we analyze the fermionic matrix models in what follows.
At low temperatures, we can insert (3.26) into (2.3) and take λ 0 to be localized at its saddle value λ 0 = 1/4γ. We arrive at an effective low temperature theory: where D φ(τ ) runs over the space of non-constant periodic functions φ(τ ). Thus we have at low energy a degree of freedom φ(τ ), with an ordinary kinetic term. Following the discussion of [25,26], this can be viewed as a low energy hydrodynamic mode (3.28) which evaluates to: Once again, we find agreement with the zero temperature correlator of the fermionic theory: where as before the fermions are anti-periodic in Euclidean time. It is convenient to consider γ > 0 such that the quartic term has the opposite sign from the vector case previously studied. At finite N , our expressions will be analytic functions of γ such that we can analytically continue γ over the complex plane. As before, we set β = 1 unless otherwise specified, in the end everything will depend on the dimensionless combination γ/β. Ai ψ iBψBj ψ jA , plus normal ordering terms which will be quadratic. The models are considerably more intricate than their vector counterpart and we will analyze them entirely in terms of the corresponding bosonic path integrals.

Effective action
Prior to integrating out the Grassmann matrix ψ iA (τ ) we have the following Euclidean action: We wish to understand the effective action of M ij (τ ) obtained upon integrating out the ψ iA (τ ). The partition function becomes: where N is a normalization constant which we will fix shortly.
Note that both the measure over M ij (τ ) and the functional determinant are invariant under the gauge symmetry: where U ij (τ ) is a τ -dependent unitary matrix. The above gauge symmetry implies that the determinant is a function of the Polyakov loop: where P denotes path ordering. However, the quadratic in M ij part of the action resembles a mass term of the gauge field and is hence not gauge invariant. Fortunately, one can take advantage of the gauge symmetry in the following sense. We can introduce a Hubbard-Stratanovich field Λ ij (τ ) such that: Note that the M ij dependent piece of the integrand is simply the generating function of correlation functions of M ij , now viewed as a U (N ) gauge field in a theory of free fermionic matrices. We can analyze this piece more carefully: Under a change of variables, one can write the Hermitean matrix M ij as: Here µ = diag(µ 1 , . . . , µ N ) is a diagonal matrix with τ -independent elements and U ij (τ ) ∈ U (N ). We can thus express (4.7) as [24]: where [DU ] is the U (N ) Haar measure. We now make the transformation Λ = UΛU † , leaving the Λ-measure invariant. The µ i integral acquires chemical potentialsλ i ≡ tr h iΛ (τ ) for each eigenvalue 4 , giving rise to a dressed partition function. All of the U dependence appears in the dτ iU † UΛ piece. This piece is independent of the constant part ofΛ ij (τ ), which we can separate out of the trΛ 2 (τ ) term. Performing theΛ ij (τ ) integral, we obtain: (4.10) The term tr UU † 2 is analogous to the dτφ(τ ) 2 in the vector case. It is decoupled from the eigenvalue integral. The normalization constant is given by: From the above expression, we can also obtain the case forγ = −γ > 0 by analytic continuation of µ i → iµ i : where now the normalization constant becomes: In this case, what was a unitary matrix U ij in (4.10) now becomes an element of the group generated by anti-Hermitean matrices.
It is of interest to note that the eigenvalue repulsion is due to a modified version of the usual Vandermonde, that involves the sin or sinh of the difference in eigenvalues.
This is characteristic of gauge theories at finite temperature [24].
It is useful to consider some simple examples, to ensure that the integrals defined above indeed give a positive definite spectrum. Consider the case with L = 1 and N = 2, i.e. 2 × 1 fermionic matrices. Here the Hilbert space is 4 = 2 2×1 dimensional.
An explicit evaluation of the integral (4.12) gives: (4.14) For L = 2 and N = 2 we find: The states add up to 16 = 2 2×2 , which is the correct dimension of the Hilbert space.
For L = 2 and N = 3 we find: The states add up to 64 = 2 3×2 , which is the correct dimension of the Hilbert space.

Summary
Thus we see that much of the structure of the vector model is carried forward to the matrix model. Instead of an ordinary integral over a single variable, we now have a partition function given by an ordinary matrix integral. Instead of a large N saddle point value for a the single variable, we will find large N eigenvalue distributions.
Finally, instead of a single low energy field φ(τ ), we now have a low energy unitary matrix U ij (τ ). The spectral information about the Hilbert space of the fermionic model is subsumed entirely into the structure of the eigenvalue integral, rather than the low energy fluctuations of U ij (τ ).
We now explore the thermal phase structure and correlations of the fermionic matrix model.

Thermal phase structure and correlations
In this section we discuss the thermal phase structure of the matrix model. We will consider the case withγ > 0. To analyze this we read from (4.10) the potential acting on the eigenvalues: As before, we fix units where β = 1 and study the partition function as a function of α ≡ L/N , which we keep fixed in the large N limit, andγ.

High and low temperature limits
We can develop a schematic idea of how the eigenvalue distribution should look. The minima of the potential acting on the eigenvalues depends on the temperature. As in the vector model, at low temperatures there are two minima. The eigenvalues will accumulate near these two minima, and will repel each other by the Vandermonde interactions. There will be many such distributions that are solutions. For instance, all N eigenvalues might be located in either minimum or one on the left and all (N − 1) others on the right and so on. Eigenvalues in one minimum can thermally tunnel to the other. As we increase the temperature, the double well profile is lost and instead there is a single minimum at the origin. Consequently the eigenvalues will be distributed around a single minimum at high temperatures.

High temperatures
In theγ → ∞ limit the Gaussian part of the potential acting on the µ i dominates.
Thus, the partition function receives a large contribution from small values of µ i and we can approximate our partition function by expanding the cosh µ i /2 in (5.1) to quadratic order. We obtain the following Gaussian matrix integral: where Q is defined in (4.13). This matrix integral has appeared in studies of Chern-Simons theory on an S 3 [28,29]. 5 The eigenvalue distribution is connected (single cut) and centered around the origin. Its explicit form is given by:

Low temperatures
At low temperatures, i.e. in theγ → 0 limit, the eigenvalue potential develops two minima around µ i = ±1/4γ. Around these, the eigenvalue potential is approximated From the result for the eigenvalue distribution (5.3) we can read the width of the distribution.
Consider the case where all eigenvalues are located around one of the minima.
We have: µ ± 1 4γ ∈ −2 cosh −1 e g/2 , 2 cosh −1 e g/2 , g = 1 2αγ . (5.5) In theγ → 0 limit we can approximate cosh −1 e g/2 ≈ log 2 + g/2. Note that for α > 2/(1 − 8γ log 2), the eigenvalue distribution has compact support away from µ = 0. For α < 2/(1 − 8γ log 2) we cannot have an eigenvalue distribution with compact support away from the origin. More generally however, due to the repulsion of eigenvalues the lowest energy configuration will favor the distribution of eigenvalues evenly among the two minima. Whether the eigenvalue distributions are disconnected will depend on α. A small α broadens the potential, enhances the effect of repulsion and connects the two distributions. For parametrically large α the eigenvalues peak sharply about µ = ±1/4γ.
In summary, the global phase structure goes as follows. At large enough temperatures there is a single cut eigenvalue distribution located near the origin. At low temperatures, the eigenvalue distribution senses two minima in the eigenvalue potential, and may be connected (for large α) or disconnected (for small α). It is interesting to note the similarity of our phase structure to the one studied recently in [12]. We will present a detailed analysis of the phase structure in future work.

Matrix Correlation functions
Finally, we would like to briefly discuss the correlation functions of ψ iA (τ ). Following the discussion for the vector case, we must invert the differential operator G −1 ij (τ, τ ) = (δ(τ − τ ) ∂ τ + δ(τ − τ )M ij (τ )). Using a parallel argument to the non-matrix case, we find: We see that the inverse operator is naturally expressed in terms of Wilson lines. Note At large N we calculate: with δτ ≡ (τ − τ ). Explicitly, we must calculate two pieces. One comes from the U ij sector: (5.9) The above piece is analogous to the contribution coming from the quadratic φ(τ ) action, as in (3.28). At large Lγ, The dominant part of the dynamics comes from the eigenvalue piece, which in the large N limit, we can express as an eigenvalue density integral: For high temperatures, the above integral can be computed numerically using the single cut eigenvalue distribution (5.3). For Lorentzian times we take δτ → iδt.
In figure 4 we show an example G(δt). We can give an analytic approximation of G(δτ ) at high enough temperatures, or equivalentlyγ large. There, the eigenvalue distribution is close to Wigner's semi-circle law, so we can approximate (5.10) by using ρ(y) ≈ (2πt) −1 4t − y 2 . Recalling that t −1 = 2α(γ − 1/8) implies that t is a small number at high temperatures. Moreover, the range of y is O( √ t) and thus to leading order in small t we find: where I n (z) is the modified Bessel function of the first kind. This function agrees well with the numerical result. For δt α (γ − 1/8), we find it decays as G(δt) ∼ δt −3/2 . This is a significantly faster decay than the vector model at high temperatures, which only experiences a sharp decay after times of order δt ∼ 4 2N γ(γ − 1/8). A similar calculation will hold for the low temperature (sub-dominant) saddle where all eigenvalues are gathered around a single minimum and t = 1/(2αγ). Higher point functions will be given by computing expectation values of products of Wilson lines.
We hope to analyze of the matrix correlators in the β = 8γ in future work.

A More general potentials
In this appendix, we express the partition function of a vector theory with general potential V (ψ I ψ I ) as a simple integral. One begins by introducing a δ-functional: δ(Λ(τ ) −ψ I (τ )ψ I (τ )) = Dλ(τ )e i dτ λ(τ )(Λ(τ )−ψ I (τ )ψ I (τ )) . (A.1) Consequently, we can express the partition function of the fermionic theory as: As in the main text, det (∂ τ + iλ(τ )) will only depend on the zero Fourier-mode of λ(τ ). Consequently, upon performing the path integral over the non-zero modes of λ(τ ), the term i dτ λ(τ )Λ(τ ) in the effective action will produce δ-functions for the Λ n with n = 0. Hence we will remain with an integral over λ 0 and Λ 0 . We can write the partition function: The bosonic path integral becomes: There is no analogue of the large N low and high temperature phases. The large N saddle point equation is: The above equation has many solutions, but the dominant saddle lives at λ 0 = 0. As we lower the temperature, we must include an increasing number of saddles to obtain a good approximation.
Correlation functions of the fermions are given in the bosonic picture by expectation values of the following non-local operator: The correlation function exhibits the Gaussian decay observed for the γ > 0 model discussed in the main text, except that it now does so for all temperatures.