Thermodynamics of a generalized graphene-motivated (2+1)-dimensional Gross-Neveu model beyond mean field within the Beth-Uhlenbeck approach

We investigate the thermodynamics at finite density of a generalized $(2 + 1)$-dimensional Gross-Neveu model of $N$ fermion species with various types of four-fermion interactions. The motivation for considering such a generalized schematic model arises from taking the Fierz-transformation of an effective Coulomb current-current interaction and certain symmetry breaking interaction terms, as considered for graphene-type models in Ref. [16]. We then apply path-integral bosonization techniques, based on the large $N$ limit, to derive the thermodynamic potential. This includes the leading order mean-field (saddle point) contribution as well as the next-order contribution of Gaussian fluctuations of exciton fields. The main focus of the paper is then the investigation of the thermodynamic properties of the resulting fermion-exciton plasma. In particular, we derive an extended Beth-Uhlenbeck form of the thermodynamic potential, discuss the Levinson theorem and the decomposition of the phase of the exciton correlation into a resonant and scattering part.


Introduction
The application of field-theoretic models with local four-fermion interactions to the investigation of fermion mass generation by spontaneous breakdown (SB) of symmetries and related phase transitions has a long history going back to the famous BCS-Bogoliubov approach to superconductivity [1,2]. Important relativistic generalizations of this approach to strong interactions considering the dynamical generation of nucleon or quark masses due to SB of a c The Author(s) 2012. Published by Oxford University Press on behalf of the Physical Society of Japan. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by-nc/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
It is interesting to note that in lower-dimensional condensed-matter systems like polymers [12] or graphene [13], the low-energy fermion excitation spectrum is effectively described by a (quasi) relativistic Dirac equation rather than by a Schrödinger equation. Such a similarity with relativistic QFT was a strong motivation to apply corresponding path-integral bosonization techniques, developed for relativistic quark models ( [4], [5], [14], [15]), to a quasirelativistic planar (D = 2 + 1) system like graphene [16]. Since the present work will continue the latter kind of investigation, let us shortly recall some main ideas of that work. There, one starts with an instantaneous four-fermion interaction obtained by a contact approximation to an instantaneous Coulomb potential using the "brane-world" or "reduced QED scenario", proposed in Ref. [17]. Next, by performing a Fierz transformation into fermion-antifermion (hole) ff -channels yields four different types of four-fermion interactions. This then leads-by means of the Hubbard-Stratonovich trick -to four types of possible exciton bound state fields E = {σ i , ϕ i }, i = 1, 2. As a result, the four-fermion interactions are now changed into Yukawa-type bilinear (Eff )-couplings. After this, the fermion path integral becomes Gaussian and can easily be carried out. The resulting Fermion determinant, rewritten as a "Tr ln" term, then determines the effective low-energy action of collective exciton fields. Note that for the application of a suitable 1/N expansion technique, it was particular useful to consider the "multilayer " case of N = 2 N degenerate fermion species (flavors) of real spin ↑ and ↓, living on N hexagonal monolayers. Thus, fermions carry here a flavor index a = (1, ...N = 2 N ). It is particularly important that in case of graphene, there exist two "valley" d.o.f. of fermions (corresponding to two "Dirac-points") which combine with two sublattice ("pseudospin") d.o.f. to form reducible four-spinors of a (quasi) relativistic Dirac theory. This just allows one to introduce chiral 4 × 4 Dirac matrices γ 5 , γ 4 , to consider a new emergent "pseudospin-valley" chiral SU (2) pv symmetry and to study both the dynamical generation of fermion masses by SB of chiral symmetry and the resulting mass spectrum of composite excitons.
The considered graphene-motivated GN 3 model of Ref. [16] includes also small symmetrybreaking contributions due to on-site repulsive lattice interactions (as well as possible phonon-mediated interactions). Since the corresponding approximation scheme does not allow a determination of the effective coupling constants from an underlying microscopic lattice theory, this generalized GN 3 model can evidently be only considered as a schematic one. Nevertheless, because of its very general structure, we find it interesting for possible applications and specializations.
It is worth mentioning that in distinction to NJL quark models of strong interactions, the fermions of the considered generalized graphene-motivated GN 3 model are not confined. Thus, contrary to the forbidden, artifact quark decay of a pion, π → qq, a real decay of a pseudoscalar exciton E → ff is possible. 2/31 Finally, let us add some remarks on the problem of renormalization. As is well-known, the GN 2 model is renormalizable and even asymptotically free [6]. On the other hand, the GN 3 model is perturbatively non-renormalizable, but becomes renormalizable in the 1/N expansion [9]. Anyway, in the condensed-matter motivated effective GN 3 model as considered in Ref. [16], there is no real problem with renormalization. Indeed, there exists a natural finite momentum cutoff Λ ∼ 1/a with a being the lattice spacing. In graphene a ≈ 2.46Å and Λ ∼ 1/a ≈ 5 keV.
The main goal of the present work is to extend the investigations of the generalized GN 3 model of Ref. [16], performed at T = 0, and to study its thermodynamics at finite temperature and fermion density (chemical potential) beyond the mean-field approximation. This requires to calculate in the large N limit the thermodynamic potential of the model including both the dominating contribution of the fermion loop and the relative order O(1/N ) correction term of exciton fluctuations around mean fields. The necessary calculations can be performed either by the diagrammatic summation of Feynman diagrams, as done for D = (3 + 1)-dimensional NJL models of colored quarks [18], or by powerful path-integral methods at finite T , µ [19] which might be further combined with 1/N expansion techniques as in Ref. [11]. In the following we prefer to use here an analogous large-N path-integral approach for calculating the grand canonical partition function and the related thermodynamic potential at (relative) order O(1/N ). Note that besides explicit symmetry breaking in coupling constants, we allow also for explicit breaking of chiral symmetry by a bare fermion mass term m 0 . This avoids dangerous IR-divergencies [11] and, most importantly, excludes the validity of the Mermin-Wagner-Coleman (MWC) no-go-theorem [20]. As is well-known, the latter forbids SB of a continuous symmetry for a (2 + 1) -dimensional system at finite temperature.
The paper is organized as follows. In Sect. 2 we shortly review the generalized, graphenemotivated GN 3 model of Ref. [16] to be further studied in this paper, including necessary definitions and motivations. Sect. 3 starts with the path-integral representation of the grandcanonical partition function and provides via the Hubbard-Stratonovich trick and after fermion integration the effective bosonized action of four types of bound-state exciton fields E = {σ i , ϕ i }, i = (1, 2). The path integral over exciton fields becomes trivial, when evaluated in the large-N (saddle point) approximation for constant "mean-field" values, and leads to the mean-field approximation of the thermodynamic potential in subsection 3.1. The respective mean-field valuesσ i ,φ i are determined from the stationarity (gap) equations as given by the minima of the thermodynamic potential. In corresponding subsubsections, we first consider the vacuum gap equations at T = 0, µ = 0 and then later the gap equation at finite temperature and chemical potential. In Subsect. 3.2 we calculate the contribution of exciton fluctuations to the thermodynamic potential at (relative) order O(1/N ) in Gaussian approximation. They are given by corresponding determinantal terms of inverse exciton propagators. The required polarization functions, entering the exciton propagators, are determined by fermion loop diagrams and calculated by standard Matsubara summation techniques. In Sect. 4 we present general expressions and equations for the mass spectrum of exciton bound and resonant states. Particular attention is drawn to the possible formation of light (would-be) Goldstone excitons. They are generally expected to dominate the thermodynamics of a fermion-exciton plasma at low temperature. On the other hand, fermions 3/31 become light due to restoration of chiral symmetry at high T and are expected to dominate in that region.
Finally, Sect. 5 contains the derivation of the celebrated Beth-Uhlenbeck (BU) formulation for the thermodynamics of a strongly correlated fermion system (comp. with [18], [21]) which we apply for a fermion-exciton plasma considered here. Note the first application of the BU formula for plasma physics using phase shifts for the two-body scattering problem of the Coulomb interaction was given by Ebeling [22]. A generalization of the BU approach for excitons in dense electron-hole plasmas was developed by Zimmermann and Stolz (see [23] and references therein) and taken over to the problem of bound state formation in nuclear matter in Ref. [24]. In the present work we will restrict ourselves to the low-density limit. Using the pole approximation, we also quote the pressure of a gas of free excitons and consider a first departure beyond such an approximation based on a Breit-Wigner type ansatz for the density of states. Moreover, similarly as in Refs. [18], [25] and [26], we indicate here on the necessity of a decomposition of the scattering phase shift into a resonant and a scattering part in order to fulfill the Levinson theorem. In Sect. 6 we summarize our main results and present the conclusions. Technical details like the algebra of Dirac matrices in D = 3, the explicit solution of the vacuum gap equations and the calculation of the exciton polarization functions are relegated to three Appendices.

Generalized GN-model: definitions and motivations
In the following we shall consider a planar system of N fermion species (flavors) described by a generalized, graphene-motivated (2 + 1)-dimensional Gross-Neveu (GN)-type of model, analogous to that considered in Refs. [8,16].
Using euclidean notations, the Lagrangian L e , L e = −L| t→−iτ , is given by a sum of two terms were L 0 is the free part including a (small) bare fermion mass m 0 and being supplemented by a chemical potential µ, and L int describes a generalized local four-fermion interaction with different coupling structures given by 1 Note that the fermion field ψ(x) in the Lagrangian L e transforms as a fundamental multiplet of the flavor group U (N ), i.e. ψ(x) = {ψ a (x)} with a = 1, ..., N . Moreover, each component ψ a (x) of this multiplet is taken to be a four-component Dirac spinor that transforms according to a reducible four-dimensional representation of the rotation group SO(3).
The reason for this choice is our interest in chiral symmetry and its possible dynamical breakdown in planar systems like graphene, which require the existence of a chiral γ 5 matrix. As is well-known, such a matrix can only be obtained in (2 + 1) -dimensions by using a reducible 4 × 4 representation of Dirac matrices [9,27]. Note further that G i , H i are general coupling constants, and γ µ (µ = 1, ..., 3), γ 4 , γ 5 , γ 45 are Dirac matrices discussed in detail in Appendix A.
Since the main purpose of this paper is to study the planar many-body system of fermions and composite excitons at finite temperature and fermion density (finite µ), we find it convenient to use from the very beginning euclidean notations for the Lagrangian, coordinates x µ = (x 1 , x 2 , x 3 = τ ) = ( x, x 3 ) (it = τ ) and γ-matrices. Finally, let us remark that spinor and U (N )-group indices are not given explicitly and summation over them is implied everywhere; moreover, the index at L e will now be omitted.
Note that the kinetic part of the free Lagrangian L 0 is invariant under discrete chiral transformations and under the discrete operations of spatial inversion P, time reversal T and charge conjugation C. Recall that in the reducible (four-dimensional) representation of the Clifford algebra, the definition of the discrete transformations P, T , C is not unique. Bearing in mind the condensed-matter roots of the graphene-motivated GN 3 model considered above, the discrete symmetry operations will suitably be taken here in such a way that they are consistent with the underlying tight-binding model on the honeycomb lattice (see Ref. [28]). By slightly modifying the corresponding definitions of Ref. [28] when using the imaginary time formalism of equilibrium thermodynamics, euclidean notations for Dirac matrices and Grassmann fields for fermions as required in path integrals (for explanations and details, see App. A.3), one has T : C : Here the spin matrix σ 2 in T acts on the indices of physical spin included in the combined flavor index a suppressed here, and the superscript t denotes transposition. (Note that in euclidean space T does not change the sign of euclidean time x 3 .) The respective transformation properties of various fermion bilinears under discrete operations P, T , C of spatial inversion, time reversal, charge conjugation and under chiral transformations with γ 5 , γ 4 are shown in table 1. It then follows that L 0 is invariant under P, T and C (with the restriction µ → −µ for C). Moreover, since L int depends quadratically on fermion bilinears, it follows from table 1 that also the total Lagrangian L is invariant.
Taking certain symmetry restrictions like It is worth remarking that in the case of a continuous chiral transformation U γ5 (1) fermion bilinears transform as follows, Analogous relations hold for continuous chiral γ 4 transformations.
Clearly, the employed approximation scheme does not allow a determination of the effective coupling constants G 1,2 arising from an underlying microscopic lattice theory, eventually including lattice vibrations. Thus, the above four-fermion interaction (14) can only be considered as a schematic one. In the following we shall therefore consider its most general form (3) and specify, for illustrations, to more symmetric cases, if necessary. Before concluding this subsection, let us yet add a remark on Lorentz-invariance. Clearly, due to the suppression of spatial components of currents and related retardation effects in the considered four-fermion interaction in graphene, the chiral invariant Coulomb expression (12) and its Fierz transformation is not Lorentz-invariant. However, this does not matter here, since (discarding again resulting vector/axial-vector interaction terms), the Fierz transformation of a corresponding Lorentz-invariant Thirring-like vector interaction reproduces, up to a factor 3, the structure (14) ( G = 0) (see, e.g., Eqs. (A.8), (A.9) of App. A in Ref. [16]).

Grand canonical partition function
The thermodynamics of the considered generalized GN-model at finite temperature T and chemical potential µ is most conveniently described by the grand canonical partition function Here H is the Hamiltonian of the system, Q =ψγ 3 ψ is the charge operator and β = 1/T with T being the temperature. It is well-known that Z(T, µ) can be expressed as a path integral, where L is the model Lagrangian (1)-(3), and denotes integration over the imaginary time (τ = it) interval τ ∈ [0, β = 1/T ] and the plane with area A = l 2 . Note that the fermionic Grassmann fields in Eq. (17) satisfy antiperiodic boundary conditions ψ( x, 0) = −ψ ( x, β).
As usual, performing the fermionic path integral (17) with four-fermion interactions requires a Hubbard-Stratonovich transformation by introducing composite scalar (σ 1 , σ 2 ) and pseudoscalar (ϕ 1 , ϕ 2 ) exciton fields. The transformed Lagrangian of the generalized GN-model then takes the form Obviously, by inserting the field equations for exciton fields (i = 1, 2) back into expression (19) reproduces the original form of the Lagrangian L in Eq. (1) 3 . Note that for a constant field σ 1 ≡ σ 1 , m = m 0 + σ 1 can obviously be interpreted as a dynamical fermion mass. It is worth remarking that taking the grand canonical expectation values (condensates . . . ) of both sides of Eq.(20) relates condensates of exciton fields (σ i , ϕ i ) to condensates of corresponding fermion bilinears ψ Γ A ψ , Γ A = Γ i , Γ i . 4 Using Eqs. (17), (19) and performing a field shift σ 1 → σ 1 − m 0 , the grand canonical partition function now takes the form (κ = m 0 /G 1 ) where bosonic exciton fields must satisfy periodic boundary conditions, The resulting Gaussian path integral over fermion fields can easily be performed and leads to a fermion determinant where the Tr-symbol means the functional trace of an operator both in euclidean coordinates and Dirac space, respectively. The grand canonical partition function is now expressed in terms of (bosonic) exciton fields alone and reads

Mean field approximation
As a first step, we now calculate the determinantal term −Tr ln S −1 for constant exciton fields, as considered in the large N (mean field -"mf") saddle point approximation of the path integral. The calculations are conveniently done by Fourier transforming the fermion fields in the respective path integral (see e.g. [29]) where the index α refers here to the (reducible) four-component Dirac spinor, ω n = (2n + 1)π/β, n = 0, ±1, ±2, ... are fermionic Matsubara frequencies and ψ αn ( p) is dimensionless. Obviously, in the mf-(large N saddle point) approximation the path integral (24), evaluated at constant mf-values (σ k , ϕ k ) mf , becomes trivial and the grand canonical partition function factorizes. The relevant thermodynamic quantity, closely related to Z, is the thermodynamic potential per area A = l 2 , Its mean-field expression can be separated into an excitonic condensate part and a fermionic contribution , respectively, where and is the inverse mean-field fermion propagator in ( p, ω n ) space. Obviously, Ω mf = O(N ) for large N . Note that the exciton mean field values (σ k , ϕ k ) mf ≡ (σ k (T, µ),φ k (T, µ)) at finite 9/31 temperature T and chemical potential µ are obtained as minima (stationary points) of the thermodynamic potential per flavor species with Γ k , Γ k being defined in Eq. (20). From Eq. (20) it follows that the mean field val- Let us next consider the fermion contribution (29). Performing the functional trace Tr which is now defined as a sum over Matsubara frequencies ω n , two-momenta p and including the trace tr over Dirac indices, i.e.
we obtain Here ε i are the four eigenvalues of the 4 × 4 matrix S −1 mf which are given by (see [8,16]) Then, Ω f takes the form where and are energies or masses associated with two (correspondingly projected) quasiparticle fields, Using standard techniques for performing the Matsubara summation of the logarithmic terms in Eq.(38) (see e.g. [29]), we arrive at the final expression . (41) 3.1.1. Vacuum gap equations. Obviously, the first term in the square bracket of Eq.(41) is just the vacuum contribution at T = 0, µ = 0, whereas the other terms describe the contributions of a free gas of two types of N quasiparticles (anti-quasiparticles) of masses M 1 and M 2 at T and µ. Note that the vacuum contribution contains an UV-divergent momentum integral which will be regularized by a cutoff Λ, using polar coordinates for momenta, i.e. | p| < Λ. The regularized vacuum contribution to the thermodynamic potential Ω vac (T = 0, µ = 0), then takes the form (comp. [8,16]) where, supposing M k Λ, a term |M k | 3 O Mk Λ was omitted. It is convenient to eliminate the cutoff parameter Λ in Eq. (43) by introducing renormalized coupling constants G r k (μ), H r k (μ), defined by corresponding normalization conditions at the normalization pointμ (k = 1, 2) Here V ≡ Ω reg vac (0, 0) denotes the effective potential, and the second derivative in Eq. (44) has to be taken atσ k =μ,σ j = 0 (j = k),φ k = 0 etc. Using Eqs. (44) and (45), one can consider new coupling-type parameters g k , h k , Obviously, the bare coupling constants G k (Λ) and H k (Λ) (considered now as functions of Λ) are independent of the normalization pointμ. It thus follows from Eqs. (44)-(47) that the new coupling-type parameters g k and h k are also independent ofμ and, using analogous arguments, also independent of Λ. For brevity, the quantities g k and h k will be in the following called coupling constants.
Eqs. (49), (50) generalize the results for the usual GN-model (see first paper of Ref. [10]) to the many-coupling case. In particular, we find that the quantities G * r k = H * r k = π/(2μ) correspond to nontrivial UV fixed points of the generalized, graphene-motivated GN-model considered here. Now, having demonstrated that the expression Ω ren vac (0, 0) is UV-finite, it is clear that the condensate valuesσ k (0, 0),φ k (0, 0) at T = µ = 0, following from the gap equations (31), (32), are described in terms of finite parameters g k , h k . Due to the explicit symmetry-breaking mass-like term κ in Eq. (48), the condensateσ 1 depends not only on g 1 , but also on κ. Finally, it is worth to remark that in condensed-matter systems like graphene, the cutoff parameter Λ is an effective quantity, restricting the low-energy region of applicability of the considered model. In particular, Λ = O(1/a) is of the order of the inverse lattice spacing a.
For illustrations, let us now consider the interesting case of an effective Coulomb interaction alone (compare discussions below Eq. (14), where now G 1 = G 2 = 1/4G C ≡ G so that g k = h k = g. Moreover, supposing G > π/Λ , the renormalized coupling is attractive, g < 0. As shown in Appendix B, in this case we getσ 2 =φ k = 0. The nonvanishing condensateσ 1 satisfies the gap equation (B9) and is given by where we introduced the mass-like term M = π/|g|. Note that the chosen positive solution realizes the total minimum of Ω vac . (η k = ±1 for k = 1, 2) Here f ± (E) are Fermi-Dirac distribution functions, and sign(x) is the sign function. For illustrations, let us again consider the particular case of symmetric coupling constants g k = h k ≡ g, admitting, however, an explicit symmetry-breaking term κ = m 0 /G related to a finite bare fermion mass m 0 . From Eqs. (53), (54) we then getφ i = 0,ρ =σ 1 in analogy to the vacuum case (comp. App. B). Let us in more detail consider the equations for the condensatesσ 1 ,σ 2 . By adding or subtracting Eqs. (53) and (55), one has (M 1,2 =σ 2 ±σ 1 ; M = π/|g| ) Comparing Eq. (58) with (57), one finds the solution −M 2 =σ 1 −σ 2 = M 1 =σ 1 +σ 2 , E p , and thusσ 2 = 0, M 1,2 = ±σ 1 . Finally, performing the momentum integral in Eq. (57) gives further Note that in the case of a small bare fermion mass m 0 (small κ), the system is characterized by a critical curve corresponding to a phase with a bare fermion mass and a phase with a larger dynamical mass m =σ 1 . In particular, by settingσ 1 = 0, one obtains for the chiral limit κ = 0 from Eq. (59) an analytical expression for the critical line of a second-order phase transition, as sketched in the phase portrait of Fig. 1 (comp. Ref. [9]), 13 As shown in Ref. [8], on the critical curve, the temperature and chemical potential are related by where in the chiral limit M = π/|g| is the dynamical fermion mass m for the vacuum T = µ = 0 (comp. Eq. (51)). Obviously, the function µ(T ) in Eq. (61) vanishes for K(T ) = 1 defining the critical temperature T c ,

Exciton fluctuations around mean fields
Up to now we have treated the grand canonical partition function (24) in the mean-field (large N saddle point) approximation. Let us now consider exciton fields as fields fluctuating 14/31 around their mean-field values, restricting us to the caseσ 1 ≡ σ 1 mf = 0,σ 2 =φ 1 =φ 2 = 0 5 . Then we have We can then decompose the inverse fermion propagator S −1 (x, y) in Eq. (22)) as where m =σ 1 . In the next step we will perform the path integration in Eq. (24) over fluctuating exciton fields in the Gaussian approximation (assuming again a suitable regularization of the vacuum term). For this aim, one has to expand the logarithm in the exponent up to the second (Gaussian) order in the fields Obviously, the linear terms in the fluctuating fields in Eq. (67) vanish either due to the gap equation (31) or due to the corresponding vanishing of Tr S mf [Γ i ] and Tr S mf [ Γ i ]. Performing then the field integration in the (bosonic) Gaussian integral leads to the following factorizing contribution from exciton correlations corresponding to σ i and ϕ i channels of interactions, Here are the inverse (bosonic) exciton propagators containing the polarization function Π E of the exciton fields E = (σ i , ϕ i ), given by corresponding fermion loops, and suitable coordinate dependence is understood. The calculation of Eq. (70) is conveniently done by using Fourier transformation to momentum-frequency space p = ( p, ω n ) and q = ( q, ν m ) for the meanfield propagator S mf ( p, iω n ) and the polarization function Π E ( q, iν m ), where ω n = (2n + 5 In Refs. [11] there were considered additional excitonic corrections to the gap expression. Such corrections were, however, shown to contribute in higher order of 1/N and will thus be omitted here. where and we found it now convenient to express the frequency dependence of the propagator by iω n . A straightforward, but lengthy calculation applying Matsubara summation techniques to loop diagrams of Fig. 2 (see App. C for details) leads to the following expression for the polarization function Π E ( q, iν m ) Here the minus sign in T ∓ refers to the σ 1,2 excitons, whereas the plus sign refers to ϕ 1,2 and f ± (E) are the Fermi-Dirac distribution functions introduced in Eq. (56). Obviously, the thermodynamic potential for exciton fluctuations follows then from the Fourier-transformed expression (69) to be given by where the Tr -symbol refers to summation in ( q, ν m ) space. Consider now the total thermodynamic potential per flavor species

. Bound states
For studying the corresponding equations for excitons in planar systems, we will closely follow here the analogous techniques of Refs. [18], [19], [25], used for studying (3 + 1)-dimensional models of colored quarks. The masses M E of exciton bound states are determined by the poles of the exciton propagators, or equivalently, from the vanishing of the inverse propagators of excitons taken at rest, q = 0, where G E are generic coupling constants (G E = {G 1,2 ; H 1,2 }), and the polarization function Π E is analytically continued to real frequencies (energies) ω, iν m → ω = M E + iη. For q = 0 one has k = p, and the only nonvanishing contribution in Eq. (73) emerges for s = s so that we can put s = −s and ss = −1 in the following. The factor T ± then becomes and the polarization function is equal to Using the relation performing the summation over the sign variable s and inserting the resulting expression into Eq. (76), one finally arrives at the (regularized) bound state equation This equation determines the exciton masses M E as functions of the coupling strength G E , cutoff Λ, temperature T , chemical potential µ and dynamical fermion mass m(T, µ). In general, it has to be solved numerically. However, it is interesting to present its solution also analytically in a closed form by using, for convenience, regularized integrals and bare coupling constants and postponing more detailed renormalization to a subsequent work devoted to the numerical analysis.
Let us start with the gap equation (31) which for the caseσ 2 =φ i = 0 takes the form Supposing a small bare fermion mass m 0 (small κ) and thus neglecting terms O(m 2 0 ), a straightforward calculation leads to the following mass spectra for the parity doublets 17/31 (σ 1 , ϕ 1 ), (σ 2 , ϕ 2 ), Here −1 is the square of the induced fermion-exciton coupling constant appearing in the numerator of the exciton propagator in the pole approximation 6 (see, e.g., Eq. (30) of Ref. [30] for the case T = µ = 0); Note that g 2 E behaves for large N as g 2 E = O(1/N ) so that the exciton masses in Eqs. (82) -(85) remain finite.
As expected, in the symmetric case of Coulomb-induced four-fermion interactions with G 1 = G 2 = H 1 = H 2 , the masses of the (σ 1 , σ 2 ) and (ϕ 1 , ϕ 2 ) excitons are degenerate. In this case we have approximate continuous chiral symmetries of γ 5 , γ 4 transformations, explicitly broken by the bare fermion mass term m 0ψ ψ. The ϕ 1 , ϕ 2 excitons are then the would-be Goldstone bosons of the explicitly broken chiral symmetry.
As discussed in Sect. 2, there exists also the alternative case of a direct symmetry breaking in coupling constants. Indeed, the inclusion of a small on-site repulsive interaction term, arising in the case of a graphene lattice, as well as of phonon-mediated interactions, leads to an additional symmetry-breaking interaction which altogether gives a modification of effective coupling constants with G 1 = G 2 = G 1 , H 1 = H 2 = G 2 , but G i = H i (see Eqs. (14), (15)). In this case the masses of the doublets (σ 1 , σ 2 ) and (ϕ 1 , ϕ 2 ) remain (approximately) degenerate, but (ϕ 1 , ϕ 2 ) excitons do not longer have typical masses of order O(κ) expected for standard would-be Goldstone bosons.
Finally, it is worth mentioning that there exists the important Mermin-Wagner-Coleman (MWC) no-go theorem [20] which for (2 + 1)-dimensional systems forbids spontaneous symmetry breaking (SSB) of a continuous symmetry for the case of finite temperature. Then, in the case of an exact chiral symmetry with κ = m 0 G 1 = 0, this requires to go beyond the leading O(1/N ) exciton contributions (see Eq. (75)) obtained in the large N expansion. In our case, this would require to take into account phase fluctuations of exciton fields related to vortex excitations and the Kosterlitz-Thouless transition [9,31]. The consideration of corresponding higher order 1/N corrections for the massless case m 0 = 0 is, however, outside the scope of this paper. In this context, it is worth to emphasize that in the case of 6 The expression for the exciton propagator contains the squared coupling constant g 2 E and presents more precisely the T −matrix expression instead of a pure propagator. Nevertheless, for brevity, it is called propagator here since both objects share the essential nonperturbative analytic property of the pole and branch cut structure in the complex plane.
18/31 explicit symmetry breaking, as considered in this paper, there do not appear dangerous infrared singularities in loop diagrams due to the lack of massless particles. Therefore, the MWC-theorem does not apply here (comp. with discussions in Refs. [11]).

Resonant exciton states
Clearly, an exciton state becomes a resonant state, if its mass exceeds two times the dynamical fermion mass so that it can decay into its constituent fermion-antifermion pair. In this case we are seeking a pole of the Fourier-transformed exciton propagators D σi , D ϕi , analytically continued into the second sheet of the complex frequency plane, iν n → z = M E − iΓ E /2. Thus, we have to solve the equation where ε E = {2m; 0} for E = {σ 1 (σ 2 ); ϕ 1 (ϕ 2 )}.
Here the two integrals I 1 , I 2 are defined by and z = M E − iΓ E /2. . Then, by separating Eq. (87) into one for the real part and another one for the imaginary part, one obtains The calculation of ReI 2 and ImI 2 appearing in Eqs. (90) and (91) can be done by using the standard formula 1 where P.V. means principal value. After some lengthy but straightforward calculations, analoguous to those in Ref. [32], one gets where the UV-cutoff Λ was used in performing the momentum integration and Θ(x) is the step function. 19/31 In deriving Eqs. (90) and (91) it was assumed that the imaginary part Γ E occuring in the argument z of I 2 (which is the full width at the half maximum of the Breit-Wigner distribution) is small compared to the mass M E and may be neglected so that the equations for M E and Γ E decouple (compare with the discussions in Ref. [18]). In the symmetric case G i = H i one expects on general grounds that with increasing temperature the constituent fermion mass m(T, µ) and thus the two-particle threshold 2m(T, µ) lowers, whereas the mass of the would-be Goldstone boson M ϕ1 (T, µ) is chirally protected and thus rather temperature independent. Therefore it is inevitable that a temperature exists where the mass of ϕ 1 coincides with that of its constituents and thus the binding energy vanishes. This temperature is called Mott temperature T Mott and is defined by the relation The effect of bound state dissolution caused by a lowering of the continuum edge in a hot and dense medium is called Mott effect due to its analogy with the physics of the insulator-metal transition according to N. Mott [33,34].

Beth-Uhlenbeck approach to the fermion-exciton plasma
Let us now in more detail investigate the thermodynamic potential Ω fl (T, µ) for exciton fluctuations, given in Eq. (74), by rewriting the logarithm of the inverse exciton propagator as an integral considered in the limit ε → 0, and Note that the frequencyindependent constant C = ln[βN/(G E ε)] compensates the logarithmic divergence of the integral at the lower boundary ε and reproduces the logarithmic dependence on N and β of the l.h.s. Next, it is convenient to rewrite the exciton propagator D λ E with coupling constant λG E as dispersion relation with the spectral function A λ E ( q, ω) Recalling that the spectral function A λ E is given by the doubled imaginary part of the analytic continuation of the propagator to the real ω-axis at positions ω ± iη, where it has branch cuts, one obtains It is worth remarking that the argument of the logarithm in the obtained expression is the very definition of the scattering matrix S E in the Jost representation [18], 20/31 Because above the threshold for elastic scattering |S E | = 1, one has where φ E is the phase of D E ( q, ω + iη). Then, using Eqs. (96)-(100), one gets where C = ln(βN/G E ). Finally, by inserting the expression (101) into the fluctuating part of the thermodynamic potential (74) and performing the trace over momenta and bosonic Matsubara frequencies ν m (compare Eq. (C11) in App. C), one obtains where g(ω) denotes the Bose-Einstein distribution function and an additive term O(ln N ) independent of m is discarded (comp. the discussion below Eq. (75)). It is convenient to split the integral (102) in the following way Next, by using the properties φ E ( q, −ω) ≡ −φ E ( q, ω) and g(−ω) = −(1 + g(ω)), and performing a partial integration, we obtain the final form of the excitonic thermodynamic potential which is the so-called extended Beth-Uhlenbeck equation. The connection with the original Beth-Uhlenbeck expression for the second virial coefficient in terms of the two-body scattering phase shift can be found in Refs. [21]. Combining a possible pole term of the S-matrix and the scattering contribution, gives Note that in the pole approximation only the first term in Eq. (106) is taken into account with E E = q 2 + M 2 E . In this case one obtains from Eq. (105) the thermodynamic potential 21/31 for a (quasi)relativistic Bose gas of non-interacting excitons Correspondingly, the pressure of a gas of free excitons is given by Here the divergent momentum integral of the excitonic vacuum energy has again to be regularized by a cutoff and should be subtracted afterwards (see discussion below). It is generally convenient to perform the variable transformation ω = q 2 + s in Eq. (105) and to introduce the density of states This leads to the alternative representation The total thermodynamic potential per flavor species is then given by Eqs. (75), (110) .
As is well-known, the equation of state and relevant thermodynamic relations are encoded in Ω(T, µ; m(T, µ)), which corresponds to the pressure density up to a sign. Obviously, only the pressure density relative to the physical vacuum term Ω phys vac at T = µ = 0 can be measured. Here the mean-field vacuum expression Ω mf vac is given in Eq. (42) and the term Ω fl T =µ=0 has to be taken from Eq. (110). Following an analogous procedure as for (3+1)-dimensional models of colored quarks in Ref. [18], the above vacuum term should be subtracted to give the physical thermodynamic potential which is again denoted by Ω (T, µ) . Clearly, such a removal of the divergent vacuum energy contributions leads to a vanishing pressure P (T, µ) = 0 at T = µ = 0. This procedure is analogous to the "no sea" approximation in relativistic mean-field approaches as, e.g., used in the Walecka model [19]. Finally, let us return to the expression given in Eq. (110) and consider the first departure beyond the pole approximation used in Eq. (107). Clearly, not too far above the Mott temperature T Mott it is justified to describe D E (s) by a Breit-Wigner type function (see a related proposal for quark matter in Ref. [19]) Here M E is the exciton pole mass, Γ E is the corresponding width and a R is a normalization factor. Below T Mott , where the spectral broadening Γ E (T ) of the state vanishes, the above expression becomes a delta function, typical for the spectral contribution of an exciton bound

Summary and conclusions
In this paper we have extended the investigation of the generalized, graphene-motivated (2 + 1)-dimensional GN-model of Ref. [16], performed at T = 0, to the case of finite temperature and fermion density (chemical potential) and have been going also beyond the mean-field approximation. The schematic model with four different types of four-fermion interactions considered here arises in a natural way by Fierz-transforming an effective electromagnetic Coulomb interaction onto four types of ff -channels. The model includes further a small symmetry-breaking repulsive on-site interaction term and a small bare fermion mass which explicitly breaks chiral symmetry. The bare fermion mass term avoids IR-divergencies and the validity of the MWC-no-go-theorem. Throughout we have used the powerful tool of large-N path integral techniques to derive the bosonized version of the grand canonical partition function and the related thermodynamic potential in terms of four composite exciton fields E = {σ i , ϕ i } (i = 1, 2). We then first employ the large-N saddle point approximation to the path integral which is based on the mean-field values {σ i ,φ i }. They are determined from the minimum of the thermodynamic potential in the form of gap equations which are studied in some detail. The mean-field approximation just leads to the dominant contribution in the thermodynamic potential and the related equation of state. On the other hand, as it is well-known from four-quark models of strong interactions (see e.g. [11], [18], [19]), light composite (would-be) Goldstone mesons (pions) play an important role in the low temperature region of a quark-meson plasma, where quarks are heavy due to a nonvanishing quark condensate. At higher temperatures the dynamical quark mass decreases, the quark threshold lowers and mesons are expected to become heavier resonances. Thus, in the high temperature region fermions should dominate thermodynamics. As emphasized in the literature, it is just this different behavior in 23/31 the low and high T regions that makes it necessary to extend the large-N path integral approach beyond mean field by including contributions of fluctuating fields at (relative) order O(1/N ) ( [11], [18]). We followed here the analogous line of argumentation for the chosen graphene-motivated GN 3 model. In particular, we present the thermodynamic potential of the resulting fermion-exciton plasma at next order in the 1/N expansion. Besides of general expressions for the fermion and exciton spectra at finite temperature and densities, one of the main results is the presentation of the exciton contribution to the thermodynamic potential in a form which extends the results of Beth and Uhlenbeck as well as Dashen et al. [21] to the in-medium case. The numerical calculation of the thermodynamic potential and the numerical estimate of the importance of the O(1/N ) corrections within the underlying generalized GN 3 model are under way.
of a symmetry group SU (2) of the free Lagrangian L 0 in Eq. (2) with group algebra where ε ijk is the Levi-Civita symbol.

A.2. Relations to graphene
For possible applications to planar systems in condensed matter physics, it is worth to compare with corresponding notations used in the chiral limit (m 0 = 0) in graphene (see [16], [28]; for an excellent review, see also [13] and references therein). As is well-known, the above SU (2) group describes there an emergent global continous "pseudospin-valley" symmetry. Moreover, one considers chiral spinors ψ ± = 1 2 (1 ± γ 5 ) ψ which have chirality eigenvalues ±1, They correspond to fermion excitations at the two inequivalent "Dirac points" K + , K − = − K + at the corners of the first Brillouin zone, where band energies vanish, ε ± ( K ± ) = 0. In particular, by performing a low-momentum expansion k = K + ( K − ) + p around these points, one obtains a linear dispersion law ε ± = ±v F | p| for massless quasiparticles on the hexagonal graphene lattice, enabling a quasirelativistic description by a Dirac Lagrangian. Definite chirality eigenvalues ±1 coincide with the so-called valley index η = ±1, characterizing the Dirac points K + and K − . Besides of two "Dirac point (valley)" degrees of freedom, graphene electrons possess in addition two "sublattice" or "pseudospin" d.o.f. . The latter ones are associated to the two existing triangular sublattices of carbon atoms in the hexagonal graphene lattice to which an electron belongs. By combining the corresponding two valley and pseudospin (sublattice) d.o.f. into a four-spinor, makes it just possible to define chiral 4 × 4 matrices (A2) and to introduce various types of fermion interactions. Finally, these analogies then allow us to apply nonperturbative methods of relativistic QFT, extensively used for studying dynamical mass generation and phase transitions for quark models of strong interactions in dependence on various types of symmetry breaking scenarios.

A.3. Discrete Symmetries
Let us next consider a convenient choice of discrete symmetry transformations P, T , C based on the condensed-matter roots of the graphene-motivated GN 3 -model considered in the text. We shall here follow the definitions of Gusynin et al. [28] which have to be slightly modified when using euclidean notations and Grassmann fields. Consider (reducible) four-spinor fields of graphene with components referring to two inequivalent Dirac points K + , K − = − K + and triangular sublattices with A and B sites 26/31 of a hexagonal honeycomb lattice given by 7 (A8) (Note that the sublattices A, B were exchanged in components at the K − point in order to get a suitable four-spinor notation.) For shortness, we omitted here the combined flavor index a including the considered numberÑ of hexagonal monolayers and the number 2 of physical spin degrees of freedom, i.e. a ≡ (i, σ) with i = 1, . . . ,Ñ and σ = (↑, ↓) so that a = 1, . . . , N = 2Ñ . Moreover, " * " denotes here the operation of Grassmann involution.
A.3.1. Spatial inversion P. By definition, P exchanges points x = ( x, x 3 ) →x = (− x, x 3 ) and momenta p = i ∂ ∂ x → −i ∂ ∂ x = − p. It is then clear that it exchanges A and B sites (sublattices) and Dirac points ( K + , K − = − K + ), respectively. A convenient choice for the spinor transformation by spatial inversion is given by or explicitely, For further application, let us also recall the well-known rule that Grassmann involution of a product of fields and a complex number c gives, for example, with c * being the complex conjugate of c. It is then, for example, straightforward to see that the current transforms as ψγ µ ψ →ψ γ µ ψ , with γ µ = (− γ, γ 3 ) , as shown in table 1.
A.3.2. Time reversal T . The antiunitary time reversal operation reverses momentum and spin, but not sublattices, and it interchanges K ± points [28]. Accordingly, in euclidean