Fermions at finite density in the path integral approach

We study relativistic fermionic systems in $3+1$ spacetime dimensions at finite chemical potential and zero temperature, from a path-integral point of view. We show how to properly account for the $i\varepsilon$ term that projects on the finite density ground state, and compute the path integral analytically for free fermions in homogeneous external backgrounds, using complex analysis techniques. As an application, we show that the ${\rm U}(1)$ symmetry is always linearly realized for free fermions at finite charge density, differently from scalars. We study various aspects of finite density QED in a homogeneous magnetic background. We compute the free energy density, non-perturbatively in the electromagnetic coupling and the external magnetic field, obtaining the finite density generalization of classic results of Euler--Heisenberg and Schwinger. We also obtain analytically the magnetic susceptibility of a relativistic Fermi gas at finite density, reproducing the de Haas--van Alphen effect. Finally, we consider a (generalized) Gross--Neveu model for $N$ interacting fermions at finite density. We compute its non-perturbative effective potential in the large-$N$ limit, and discuss the fate of the ${\rm U}(1)$ vector and $\mathbb{Z}_2^A$ axial symmetries.

B Fermi gas in a magnetic field: weak field expansion 29 1 Introduction It has been long appreciated that systems of free scalars and free spin-1/2 particles at low temperature and nonzero charge density have strikingly different properties.The former give rise to the phenomenon of Bose-Einstein condensation and are characterized by the spontaneous breaking of an internal symmetry.The latter fill out a Fermi sphere in momentum space, and have unbroken internal symmetry.These properties can be traced back to the statistical properties of the underlying fields, bosonic and fermionic respectively, as dictated by the spin-statistics theorem in relativistic quantum field theory (QFT).
In the presence of interactions, the fate of the internal U(1) symmetry in a state of finite charge density is less clear.For systems of interacting scalar fields, in [1] we have provided evidence that the ground state at finite chemical potential cannot develop a charge density unless it spontaneously breaks the U(1) symmetry (see also [2,3] for previous related discussions).We have proven this statement in perturbation theory at one loop for generic non-derivative scalar's self-interactions, and in a O(N ) vector model at large N .A complete proof is still missing, but these results are suggestive of the generality of the phenomenon.This connection plays a crucial role also in the superfluid effective theory approach to the large charge sector of Conformal Field Theories [4][5][6][7] (see also [8] for a review).An important alternative phase is represented by Fermi liquids [9][10][11] and one would like to understand under what conditions this phase arises.
A difficulty that one encounters when trying to address this question in full generality is that of properly accounting for the statistics of the fields [1].This should enter nontrivially to encode the different properties of bosons and fermions at nonzero charge density.In order to make progress in this direction, it seems instructive to reconsider the question about spontaneous symmetry breaking in fermionic QFTs at finite chemical potential and analyze the differences with respect to bosonic systems.In this spirit, our goal in this work is to revisit for fermionic fields some of the aspects that we addressed in [1] for scalar theories.In particular, we will pose the following question: given a fermionic theory at finite charge density, is the U(1) V symmetry realized linearly or spontaneously broken?To answer the question, we will adopt a path-integral approach and show explicitly that, to correctly compute physical quantities such as the free energy density of a relativistic Fermi gas, one needs to carefully take into account the iε term in the path integral-which allows to project on the correct ground state of the Hamiltonian at finite µ.As we will discuss, the iε term is a key ingredient that is responsible for the different analytic structure in the path integral and the resulting partition functions of fermions and bosons.Analyses of systems at finite chemical potential (or density) and zero temperature are often formulated as the T → 0 limit of finite temperature QFT calculations.The subtlelties associated to the zero-temperature limit have been analyzed in an interesting recent work [12].Motivations include the study of the phase diagram of QCD at finite density [13][14][15].In this article we take a complementary approach, and discuss how to consistently perform finite density calculations in fermionic QFTs directly at T = 0.With this understanding, we will compute the partition function of a free Fermi gas in the presence of a source term in the path integral-which will play the role of an order parameter for the U(1) V symmetry-and show that the vacuum expectation value of the order parameter vanishes in the limit in which the source goes to zero, proving that the U(1) V symmetry is unbroken for free fermions.We will then extend the simple case of the free Fermi gas in two different directions.
First, we will include a non-trivial magnetic background for the Fermi gas and compute explicitly the fermionic functional determinant in the presence of a magnetic field.Our results can have applications in the study of QED in intense magnetic backgrounds and at finite density.The dynamics of QED with strong external fields is an old subject [16][17][18], see e.g.[19] for a recent review.Applications range from the astrophysics of strongly magnetized pulsars [20][21][22] to laser and plasma physics [23,24].The quantum effective action of QED at finite temperature and density has been computed in Refs.[25][26][27].Similar results have been obtained in 2+1 dimensions [28].These works considered systems at finite temperature and express the results as thermal integrals.With our approach, working directly at T = 0, we obtain closed form analytic results, non-perturbatively in the electromagnetic coupling and magnetic field.As a warm up, at zero chemical potential, we provide a modern derivation of the Euler-Heisenberg quantum effective action in its zeta function representation [29][30][31].We then obtain closed form analytic expressions for the finite density contributions, including the de Haas-van Alphen effect.
Second, we will include interactions among fermions.We will do so by considering a 3+1 dimensional generalization of the Gross-Neveu model, for a system of N Dirac fermions transforming in the fundamental representation of a vectorial U(N ) global symmetry group.In the limit of large N , the model is solvable analytically and allows one to compute the fermionic path integral exactly.Although we will not prove in full generality that the U(1) V symmetry is unbroken at finite density, we will provide evidence that the system can support a finite density phase with unbroken symmetry.This should be contrasted with the result discussed in [1] for the O(N ) vector model, where we showed that, for a system of N scalars, the breaking of the U(1) symmetry is instead inevitable at nonzero charge density, in the large-N limit.
The work is organized as follows.In section 2 we start with some general considerations about the spontaneous breaking of the U(1) V symmetry in fermionic theories at finite density from a path-integral point of view.In section 3 we discuss the role of the iε term in the fermionic path integral, showing the difference with respect to the bosonic case, and provide a proof that the U(1) V symmetry is unbroken for free fermions.In section 4, we discuss various aspects of a Fermi gas in a magnetic background, including a derivation of the QED effective action and magnetic susceptibility at finite density.The (generalized) Gross-Neveu model is then discussed in section 5. Some technical aspects and useful relations are collected in the appendices.

Conventions:
We adopt natural units ℏ = c = 1, and the mostly minus convention for the Minkowski spacetime metric.We assume that the chemical potential µ is non-negative (µ ≥ 0), for ease of notation.The free energy density as a function of µ, is defined as To simplify the notation we will suppress factors of volume.When needed, they can be reintroduced by dimensional analysis.Conventions on γ matrices are detailed in appendix A. The symbol μ is reserved for the MS renormalization scale, and its occurrence should be clear by context.
2 Symmetry breaking at finite density?
In order to address the question of whether the U(1) V symmetry is linearly realized or spontaneously broken in a fermionic theory at finite density, we look for an order parameter for the symmetry, that is: an operator O which has a non-vanishing commutator with the charge operator If such an operator has non-zero expectation value on the ground state, i.e. ⟨δO⟩ ̸ = 0, the symmetry is spontaneously broken (i.e. the ground state is not invariant under Q).
Let us start by considering a system of free massive Dirac fermions in 3 + 1 dimensions at finite density.The underlying dynamics is controlled by the Lagrangian This theory has a continuous U(1) symmetry 3) The Noether current associated to the symmetry is corresponding to the number of left-handed plus right-handed particles (the vectorial current).
We would like to show that all the (local) order parameters quadratic in the field variables have zero expectation value, suggesting that the symmetry is linearly realized.The argument we present here is naive, since it assumes that the partition function is differentiable with respect to the Majorana sources j and j * around (j, j * ) = (0, 0).We will justify this assumption later in explicit examples of fermionic theories.
A prototypical order parameter for the vectorial U(1) V symmetry is the Majorana term In order to compute its vacuum expectation value we use the path integral formulation for fermions at finite µ and include a source term.The finite µ action will be rederived later, but for now let us consider directly the finite µ path integral with Majorana source The vacuum expectation value of the Majorana order parameter is given by Since the action is quadratic in the field variables, the path integral can be performed exactly and expressed as a determinant.The finite µ Lagrangian, including source terms, can be written as follows (up to boundary terms): where The path integral can then be formally expressed as: . (2.9) We can rewrite the determinant as follows: where we used that the spin 1/2 Dirac representation is even dimensional.We arrive at the formula: By induction on n it is easy to show that the matrix where p i are formal polynomials which can also depend on the differential operator A. Therefore, evaluating the trace, the result is of the form log Z [j; µ] = const + O (j * j) .(2.13)This implies that the vacuum expectation value in equation (2.7) is zero: In this derivation we assumed that the expansion (2.13) is well defined, or in other words that the function is differentiable in (j, j * ) = (0, 0).For the free fermion we will explicitly evaluate the first two terms of this expansion in section (3.4), verifying the validity of this assumption.
The argument we presented relies on the fact that the Majorana source term has elements only on the principal diagonal of the quadratic form in equation (2.8), whereas the kinetic term and the term describing the chemical potential are off-diagonal.It can be easily generalized to order parameters of the form δO = c 1 Ψ c Ψ + c 2 Ψ c γ 5 Ψ + h.c. .

Free fermions at finite density
To gain better control of these formal manipulations we would like to explicitly compute the partition function of a system of free fermions at finite µ and zero temperature in the path integral approach.In the case of free bosons at finite µ, the path integral with the iε prescription is convergent only for µ < m [1].It is convenient to define D µ to be a formal µ-dependent covariant derivative.In the scalar case the computation reduces to the evaluation of the (inverse) determinant of the operator −D µ D µ − m 2 , which can be easily seen to be independent of µ for µ < m.On the other hand, for µ > m the gaussian path integral is divergent and one is forced to include interactions in order to stabilize the system for values of µ > m.The situation is however different in the case of fermions: the quadratic fermionic path integral is always equal (by definition) to a functional determinant, and for a (constant and homogeneous) chemical potential one obtains again the operator −D µ D µ − m 2 , both for µ < m and for µ > m.This determinant should be independent of µ for µ < m, essentially for the same reasons as in the scalar case, but should depend explicitly on µ in the regime µ > m, corresponding to the free Fermi gas phase.How can these two conditions be mutually consistent?How does the µ dependence in the determinant arise?In order to answer these questions we shall carefully analyze the iε term that selects the finite density ground state, and compute the relativistic Fermi gas free energy in the path integral approach.
We start from the theory of a free Dirac fermion of mass m, as in eq.(2.2).We are interested in studying this system at zero temperature, in the presence of a chemical potential for the vectorial charge (i.e., at finite charge density).To describe this system we first switch to the canonical formalism and include a term describing the system at finite chemical potential.In canonical quantization the conjugate variables are The Legendre transform of the Lagrangian (2.2) gives the canonical Hamiltonian density for a free Dirac fermion: (3.2) The system at finite chemical potential can then be described by the effective Hamiltonian where Q = d 3 x Ψγ 0 Ψ.Going back to the Lagrangian formalism we obtain

iε term for the fermions
To describe a finite density phase we want to project on the ground state of the operator H µ .
In order to do so we introduce an iε term in the Hamiltonian path integral for the fermions.We add to the action functional a term iεH µ , and obtain: This is equivalent to performing a Wick rotation t → (1 − iε)t, as evidenced by rewriting the functional integral as The quadratic path integral for a complex fermion is equal to the determinant of the kinetic operator.Up to an overall constant factor, we therefore obtain both for µ > m and for µ < m.Defining for convenience a covariant derivative D α = ∂ α − iµ δ α0 , and working at linear order in ε, we arrive at

Quantum mechanical fermionic oscillator
As a warm-up, let us first consider the quantum mechanical fermionic oscillator, i.e.QFT in 0 + 1 dimensions. 1 In this case there is no γ matrix and we simply need to compute Let us consider the charge density operator Q. Taking a Fourier transform and working at linear order in ε we have Using the distributional identity and integrating we obtain the charge density The −1/2 term is an additive constant that can be renormalized away.In quantum mechanics the charge operator is proportional to the Hamiltonian (but dimensionless), so that it corresponds to nothing else but the zero point energy of a fermionic quantum oscillator: The eigenvalue of the number operator N is 1 for µ > m, corresponding to the full "Fermi sea" of a nonzero charge state, and 0 for µ < m, corresponding to the zero charge ground state.

The Fermi gas in 3+1 dimensions
In the case of a Dirac fermion in 3+1 dimensions we need to deal with the spinor structure.
Using charge conjugation and the transformation property of the γ matrices Cγ µ C −1 = −(γ µ ) T , we have the identity where we used (det C)(det C −1 ) = 1, and that the determinant is preserved by transposition.Taking the geometric average, using the fact that the Dirac matrices are even dimensional, and neglecting terms of order ε 2 we obtain: For a constant and homogeneous chemical potential we can use the identity / D / D = D µ D µ .Taking a Fourier transform and using (log det = Tr log), we arrive at where and the factor of 4 comes from the trace in Dirac space.Notice that in relating the determinant of the Dirac operator to the determinant of an operator which is a multiple of the identity in Dirac space we have already made use of charge conjugation.The resulting singularity structure is asymmetric under µ → −µ, but the conjugate contribution is already accounted for by the overall factor of 2.
Let us analyze in more detail the singularities of the integrand in the complex p 0 plane, for fixed spatial momentum ⃗ p: there is a branch cut joining the two endpoints defined by the condition It is technically convenient to slightly deform the iε term in such a way that the singularities are all shifted by the same amount, with the direction of the shift determined by (3.18).After this deformation, we have It is straightforward to check that the integrand in eq.(3.16) is single-valued once this branch cut is fixed.
Figure 1: Analytic structure and contour deformation for µ > m and |⃗ p| < p F .In the first step, we deform the integration contour from the real axis (blue) to the imaginary axis (red), see panel a).The deformation is smooth and the integrals along the two arcs in the first and third quadrants cancel each other.In the second step, we deform the contour with a shift by −µ.The integral along the vertical green line reproduces the µ = 0 path integral, while the contributions from the discontinuity across the branch cut and the segments at infinity reproduce the finite µ free energy density of the relativistic Fermi gas, see eq. (3.25).
For µ < m the branch cut endpoints p 0,± lie in opposite quadrants (the second and the fourth).Let us choose the cut so that p 0,± are joined through the point at infinity, without touching the first and third quadrants and the vertical strip comprised between p 0,± .With this choice we are free to deform the p 0 integration contour as follows: first Wick rotate from the real axis to the vertical ip 0 axis and then translate to the left by µ.Doing this we do not encounter any singularity and it is easy to check that the contributions from the contours at infinity cancel each other.We are left with an integral on the vertical line, with two singularities lying at a distance ± ⃗ p 2 + m 2 , on the left and the right.This can be recognized as the Euclidean path integral of the µ = 0 theory.Therefore we have This is reminiscent of the so-called "Silver Blaze" property of QCD at (small) nonzero isospin chemical potential [32].Notice, however, that our chemical potential is associated to the U(1) vectorial symmetry and is thus more similar to a baryon chemical potential.
For µ > m, the location of the endpoints depends on the relative size of |⃗ p| and the Fermi momentum p F : They lie in opposite quadrants (the second and the fourth) only when |⃗ p| > p F .In this case the previous analysis is unchanged.For |⃗ p| < p F , instead, both endpoints lie in the second quadrant and more care is needed.In this case we choose the cut to be the horizontal segment joining the two points, see fig. 1. Trying to follow the same contour deformation procedure as before we obtain some additional contributions due to the singularities.The Wick rotation goes smoothly and one can easily check that also in this case the contributions from the large arcs cancel each other (left panel of fig.1).The translation to the left by µ gives instead two contributions: one (I 1 ) from the discontinuity along the branch cut of the log; the other (I 2 ) from the integrals along the segments (Re(p 0 ), Im(p 0 )) = ([−µ, 0], ±iΛ), see the right panel of fig. 1.We have: where In the second equation we paid particular attention to the definition of the phase with respect to the cut.We thus arrive at: which can be recognized as the zero temperature limit of the free energy density of a relativistic Fermi gas.The integral can be computed explicitly and we obtain: (3.25)The charge density is simply given by The same result for the charge density can be obtained by first performing the derivative and then carrying out the contour integration (in which case only a pole is present, instead of a branch cut).These equations match the well-known results for the relativistic Fermi gas, usually derived in the canonical formalism.

Absence of symmetry breaking for free fermions
In this section we try to compute explicitly the partition function in presence of a Majorana source term, using equations (2.9) and (2.10) as a starting point: where The operator A is invertible, therefore we can consider the decomposition: from which it follows that (3.30) The first factor corresponds to the usual j = 0 fermionic determinant (see e.g. the previous section), and we can thus focus on the second term that will give rise to the finite j corrections.Since the operator A is independent of j, j * and the operator B is proportional to j, we see that the partition function can depend on j and j * only through the combination jj * , in agreement with the argument of section 2.
The operators B, B † can be rewritten as allowing to express the second determinant as In order to show that no subtlelties arise in equation (2.14) we want to show that the first order term in the jj * expansion of the logarithm of the partition function, log Z µ [j, j * ], is regular up to the usual UV renormalization.In other words, defining we would like to show that δ jj * [µ] − δ jj * [0] is finite.After straightforward manipulation we obtain Carrying out the trace in Dirac space and using the transpose of the γ matrices (see appendix A) we arrive at (3.34) We can now compute δ jj * [µ] − δ jj * [0] using the contour integration approach previously described, obtaining (3.35) The second term in parenthesis has a double pole in p 0,+ and has therefore vanishing residue.Evaluating the residue of the first term, we arrive at which is manifestly finite and regular.We can thus conclude that where the first two coefficients are finite.We can thus conclude that the finite density phase does not induce an expectation value for the Majorana operator, and therefore that the U(1) V global symmetry is unbroken.This explicit computation puts on firmer grounds the naive argument of section 2 in the case of free fermions.

The Fermi gas in a magnetic background
Having understood how to deal with free fermions at finite density in the path integral approach, we would like now to consider the dynamics of electrons in a magnetic background, i.e. finite density QED with an external magnetic field.

The functional determinant
Turning on gauge interactions, we consider the fermionic path integral in a magnetic background.The covariant derivative is now given by where e is the QED coupling constant and A is a background gauge field associated to a homogeneous and constant magnetic background.Following the same steps as in the previous case, the fermion one-loop contribution to the partition function is given by: where now / D / D = D µ D µ − e 2 σ µν F µν .The differential operator D µ D µ acts on a test function f (x) as: For a homogeneous and constant external field we can choose a gauge in which A 0 = 0 and ∂ µ A µ = 0, so that this expression is simplified.Let us consider the case of a magnetic field ⃗ B = B ẑ along the z-axis and choose ⃗ A = Bx ŷ.Using the explicit expression and reorganizing the terms, the partition function (4.2) is given by the determinant of the differential operator In order to find its eigenvalues we use the method of separation of variables and look for eigenfunctions of the form f (x, y, z, t) = e −ip 0 t e ipyy e ipzz g(x).
This is an eigenfunction of L provided that g(x) satisfies: with regular solutions where H n (ξ) are the Hermite polynomials.We assumed implicitly eB > 0 for notational simplicity, otherwise eB should be replaced by |eB|.We have thus found the eigenvalues of L. Denoting the spin with σ = ±1/2: Taking into account the multiplicities and the normalization of the eigenfunctions g n , the partition function is given by log Z (1) (4.10)

A modern derivation of the nonperturbative Euler-Heisenberg effective action
Before discussing the finite µ case, let us first rederive the one-loop QED effective action in a background homogeneous magnetic field.We start from the expression of the partition function (4.10), evaluated at µ = 0.In this case we can safely perform a Wick rotation to Euclidean momenta and use the identity log( The resulting expression involves an integral in d 2 p which we evaluate in dimensional regularization to obtain log Z (1) [A; 0] = − ∂ ∂α where ζ(s, q) is the Hurwitz zeta function defined by (4.12) Evaluating the derivative at α = 0, we arrive at log Z (1) From the definition of the Hurwitz zeta function we have the identity: which lets us rewrite the effective action as log Z (1) [A; 0] = −i eB 2π 2eB 4π where we denoted the (arbitrary) renormalization scale as μ, not to be confused with the chemical potential µ, and ζ ′ (s, q) = dζ(s, q)/ds.Using the identity where B n (x) are Bernoulli polynomials, the first (divergent) term can be rewritten as showing explicitly that all the divergencies can be removed by a cosmological constant counter-term and a wave-function renormalization for the magnetic field B or, equivalently, the field strength F µν , which for vanishing electric field satisfies the relation B 2 = −F µν F µν /2.We choose to work in the MS renormalization scheme and remove all the terms in (4.18).For simplicity we also chose the renormalization scale to be μ = m.
We have obtained the full partition function in closed form, non-perturbatively in the coupling e and in the non-dynamical external magnetic field B. It is convenient to express the result in terms of β = 2eB/m 2 : log Z (1) to be compared with the well-known result of Euler and Heisenberg [16][17][18] log Z (1) An equivalent result has been obtained long ago through zeta function regularization methods, see for instance chapter 6 of Ref. [29], and Refs.[30,31] for a detailed discussion.We can check explicitly that the zeta function representation agrees with the Euler-Heisenberg representation in the weak field limit, where β = 2eB/m 2 → 0. Using the asymptotic expansion (see section 25.11 of [33]): together with equation (4.17), we recover the expansion −i log Z (1) [A; 0] where F 2 = F µν F µν = −2B 2 and α = e 2 /4π, which reproduces the F µν F µν terms of the Euler-Heisenberg Lagrangian with the correct numerical coefficients.
The zeta function representation we derived in equation (4.19), allows to easily obtain the strong field limit β = 2eB/m 2 → ∞: where we chose again μ = m.We see that in the strong field limit the only dependence on the fermion mass appears in the logarithmic running term, where m plays the role of an IR regulator.The purely electric case is obtained from the substitution B → iE.The (nonperturbative) Schwinger effect [18] is encoded in the imaginary part of the full result (4.19) and is not captured by the asymptotic expansion, which is real.
Adding the tree-level contribution, the free energy for a constant and homogeneous magnetic field minimally coupled to a charged Dirac fermion is given by which (for μ = m) is with β = 2eB/m 2 .

Finite density QED in a magnetic background
Reintroducing the chemical potential we consider equation (4.10): The log in the integral has branching points at The two points p 0,± lie in the same (second) quadrant when (4.28) We can now perform the integral by deforming the contour as detailed in the free Fermi gas and write the partition function as the sum of a term given by the partition function at µ = 0 and a term that we know analytically: log Z (1) [A; µ] = log Z (1) [A; 0] + Shifting the index n by a unit for σ = − 1 2 we can rewrite the result as log Z (1) [A; µ] = log Z (1) [A where and the sum over ℓ gives a sum over Landau levels.It is convenient to express the result in terms of the free energy density The integrals can be computed in closed form where we defined (4.34)It is possible to rewrite this result in terms of logarithms using the identity arctanh In the strong field case, when 2eB > (µ 2 −m 2 ) and l max = 0, only the first term contributes and the full result takes a simple form where in the first line we assumed also 2eB ≫ m 2 .On the other hand, in the weak field limit 2eB ≪ (µ 2 − m 2 ) the sum can be recognized as a Riemann sum in the variable It is convenient to add and subtract the ℓ = 0 term, so that in the limit 2eB → 0 the sum from ℓ = 0 to ℓ max converges to the definite integral in dq 2 with extrema

.38)
The leading order result becomes more transparent by taking a step back and writing the terms of the sum as integrals in dp z , as in equation (4.30), (4.39) One can now recognize the integral as an integral in cylindrical coordinates over a domain which is a ball B µ of radius p F = µ 2 − m 2 (the Fermi momentum).We thus recover the B = 0 result of equation (3.24), up to a µ-independent cosmological constant term coming from f (B = 0, 0) that we neglect (see eq. (4.22)).The lowest order correction can be computed explicitly analytically, see appendix B for its derivation: (4.40)

Magnetic susceptibility and the de Haas-van Alphen effect
The free energy of the system is not a directly observable quantity.More interesting quantities are the magnetization of the system, or the magnetic susceptibility, describing the response of the system to changes in the background magnetic field.The latter is defined as and describes the induced magnetization response. 3Using eq. ( 4.25) we can obtain the magnetic susceptibility of the QED vacuum in the presence of a magnetic background, and an analytic derivation of the de Haas-van Alphen effect.Since χ B is a dimensionless quantity, the result takes a particularly nice form once written in terms of dimensionless variables.The natural variables are β, measuring the magnetic field strength, and the chemical potential relative to the fermion mass μ = µ/m.Using the magnetic susceptibility at finite µ can be expressed as where α is the fine structure constant α = e 2 /4π.The function g(β, μ) is a rescaled version of the quantum contribution to the free energy density (4.33), given by4   where (4.45) The last term can be rewritten using the Heaviside theta function as and one might worry of picking up δ (ℓ max − ℓ) (and δ ′ ) contributions from the derivatives in eq.(4.43).These contributions, however, vanish identically up to second derivatives in β, since G ℓ and ∂ β G ℓ vanish identically when ℓ max crosses an integer value and ℓ = ℓ max .Noticing also that G 0 depends linearly on β, we arrive at The first term describes the magnetic response of the QED vacuum at finite density, for arbitrary values of the external magnetic field.We display the zero density contribution in figure 2, in units of α/π.The weak field regime is governed by the asymptotic expansion of the magnetic Euler-Heisenberg Lagrangian (4.22).This provides a good approximation only for β < 1.Notice that, being an asymptotic expansion, the region of β for which the approximation works well shrinks to zero as we increase the number of terms.In particular, the expansion can never work well beyond β ∼ 1.In the opposite limit of strong field, from eq. ( 4.23) we obtain instead The strong field limit provides an estimate accurate at the 10% level for β ≳ 20, and at the 1% level for β ≳ 100.We see, therefore, that in order to describe the magnetic response of zero density QED in the range 1 ≲ β ≲ 20 the full non-perturbative result is needed.Consider now the second term in eq. ( 4.47), describing the finite density contribution to χ B .In the limit of very strong magnetic field limit one has ℓ max = 0, so that the magnetic response at finite density is (exactly) the same as that at zero density.For a strong field, but weak enough that at least one Landau level is filled, the finite density contribution dominates.From the expression in eq. ( 4.45) we have so that the susceptibility diverges whenever μ2 −1 is an integer multiple of β, and ℓ max jumps by one unit.This oscillatory feature, displayed in figure 3, occurs every time a Landau level achieves integer filling.This periodic property of the magnetic susceptibility χ B as a function of β is known as the de Haas-van Alphen effect [34,35].The period is given by

.50)
where ) is the area of the Fermi surface.In the weak field limit, on the other hand, using eq.( 4.40) we obtain the continuous approximation for the susceptibility As we show in figure 4, the weak field approximation is obtained only as an average.

The four dimensional Gross-Neveu model at finite density
Having understood the case of free fermions in the path integral approach, we would like to include interaction terms among the fermions and see if our conclusions are modified by the presence of interactions.To this end, we shall consider a generalization of the Gross-Neveu model [36] to d = 3 + 1 dimensions, and work in the 1/N expansion.In our analysis we make the assumption that spatial translational invariance is unbroken by the ground state (i.e., that the system is in a homogeneous phase).Spatially inhomogeneous phases are known to occur in the Gross-Neveu model in 1+1 dimensions [38] and in its chiral version [39], however there is some evidence [40] that homogeneous phases are stable in 2 + 1 dimensions.We shall assume this to be the case also in 3 + 1 dimensions.
In 3 + 1 dimensions, the Gross-Neveu model can be defined for a system of N Dirac fermions transforming in the fundamental representation of a (vectorial) U(N ) global symmetry group.We consider the massless case for simplicity, but our discussion can be extended to include a mass term.The interaction term is a four-fermion interaction: where the sum over the flavor index a is implicit.This action has a U(N ) symmetry of which U(1) V is a subgroup.The axial symmetry U(1) A is explicitly broken by the interaction term, however a discrete axial Z A 2 survives, acting on the fermions as: This symmetry forbids the mass term, and prevents its appearance under renormalization group flow.

Zero density
We can rewrite the Lagrangian as a quadratic fermionic Lagrangian by introducing a scalar auxiliary field [36].To do so we simply add a Gaussian term to the Lagrangian: which has the only effect of changing the normalization of the path integral by a constant factor.The four-fermion interaction terms are canceled by the auxiliary term, and we are left with The discrete axial symmetry Z A 2 acts on σ as σ −→ −σ. (5.5) In the large N limit, N → ∞ with λ fixed, the auxiliary field σ becomes infinitely heavy, and the σ propagator is suppressed by a factor 1/N .At leading order in the 1/N expansion we can therefore treat the scalar field σ at tree level, and consider fermionic loops only.The Gross-Neveu model is renormalizable in two dimensions, and no additional counterterms are needed beyond the four-fermion interaction, or equivalently the quadratic term for the auxiliary field σ.In four dimensions, on the other hand, the theory is nonrenormalizable, and an infinite number of counterterms would be needed to renormalize the theory to all orders.The situation simplifies however at leading order in large N : a single counterterm gσ 4 is needed to renormalize the theory, to all orders in the coupling λ.
To see this, consider now Z[σ; 0], which can be computed by standard methods, analytically continuing the momenta to Euclidean space.Up to an additive constant, we have: where p is in Euclidean space, i.e. we replaced p 0 → ip 0 .We regulate the theory adopting dimensional regularization.The divergent contribution for d → 4 is We can renormalize the theory by including of a local g 4 σ 4 counterterm in the action.In the MS scheme, neglecting the cosmological constant term, we obtain where μ denotes the renormalization scale and ḡ4 is the MS running coupling.The effective potential, as it stands, has an unstable direction for large values of σ, due to the negative logarithmic contribution enhanced by the large N factor.We are therefore forced to take into account higher dimensional operators, that stabilize the theory.For this purpose we include higher order terms in the auxiliary scalar field formulation, taking the Lagrangian as the definition of the fermionic model.Only even powers of σ appear, to preserve the discrete axial symmetry Z A 2 .For the purpose of illustration, we consider the case g 6 > 0 and g 2n = 0 for n > 3, but the results we discuss can be generalized straightforwardly.In the purely fermionic formulation this corresponds to adding an infinite series of local interactions of the form ( Ψa Ψ a ) n , with coefficients determined by λ, g 4 , g 6 and N .
At leading order in 1/N , the coupling λ does not run and the large N theory can be defined by specifying the (arbitrary) value of λ.We can then express the running quartic coupling ḡ4 (μ) in terms of a dimensionless physical coupling κ 4 , evaluated at the scale Computing the beta function of ḡ4 (μ), and solving its RG evolution, we find (5.11) Similarly, the coupling g 6 can be expressed in terms of a dimensionless coupling κ 6 by defining g 6 = κ 6 M 2 . (5.12) We shall assume that M is the only dimensionful scale in the problem, so that κ 4 , κ 6 are of order one.It is convenient to adopt units M = 1 (factors of M can always be reintroduced by dimensional analysis), and set κ 2 ≡ sign(λ).
Expressing the effective potential in terms of physical couplings we arrive at (5.13) In taking the large N limit we keep the couplings κ 4 , κ 6 fixed.Loops involving κ 4 and κ 6 vertices, and σ as an internal line, are therefore subleading in the 1/N expansion.We notice that the μ dependence of the effective potential drops out once we express the result in terms of physical couplings.This property is consistent with the fact that the minimum of the effective potential is related to observable quantities.The effective potential (5.13) at large N has a global minimum σ 0 at + O (log log N ) . (5.14) A more accurate estimate can be obtained in terms of the negative branch of the Lambert W function (neglecting κ 2 , κ 4 , which are suppressed at large values of σ 0 ) 24π 2 e κ 6 N .(5.15)This minimum always exists for λ > 0 (κ 2 = +1) and κ 4 , κ 6 ≈ 1 , whereas it requires large enough N for λ < 0 (κ 2 = −1).The large N scaling of σ 0 would be modified by the presence of higher order interactions, e.g.κ 8 ≈ 1, however the qualitative features of the model would be unchanged.The expectation value of σ breaks spontaneously the discrete Z A 2 symmetry, generating an effective mass term m eff = σ 0 for the fermions, similar to the mass generation mechanism for quarks and leptons by the Higgs field in the Standard Model.The vectorial symmetry U(N ), and in particular U(1) V , is instead unbroken.

Finite density
We can introduce a chemical potential for the U(1) V ⊂ U(N ) symmetry, in such a way that the SU(N ) subgroup is preserved: (5.16) We would like to show that this theory can support a finite density phase with unbroken U(1) V symmetry.A necessary condition for this is that the ground state energy density of the system has a non-trivial µ dependence, for µ larger than some critical value µ crit .We will show that in the large N fermionic model under consideration this expectation can be met without the need for U(1) V breaking order parameters.This behavior should be contrasted with that observed in the case of interacting scalar fields, and in particular in the O(N ) model with quartic interactions, analyzed in [1], where it is shown that the system can support a finite density phase only in the presence of a non-zero expectation value for an U(1) V breaking order parameter.We shall find that the critical value µ crit to support finite density coincides with σ 0 , the effective (pole) mass of the lowest lying charged states, as expected on physical grounds and argued also in the scalar case [1].
In the background of σ, the fermionic part of the action at finite µ is simply given by (5.17) This amounts to N copies of the free fermion (3.4), with m → σ.The finite µ fermionic path integral can therefore be computed exactly using the contour integration method previously described.We are interested in computing .18)From this we will compute the free energy density of the system by setting the auxiliary field at its background value: We already computed µ effective potential V eff (σ; 0) in eq.(5.13).All we need to compute is, therefore, the difference V eff (σ; µ) − V eff (σ; 0).Given the fermionic action (5.17), the contribution V eff (σ; µ) − V eff (σ; 0) can be computed following the same complex contour integration technique we used for free fermions in the path integral formalism, as in sec.3.3.
The result is simply given by N times (3.25), after the substitution m → σ: It follows immediately that as long as µ < σ 0 the effective potential is not modified in a neighborhood of σ 0 , the expectation value of σ at µ = 0.One can easily check that for µ < σ 0 the global minimum of the potential V eff (σ; µ) is still given by σ 0 .The fermions have zero charge density, since and For µ ≥ σ 0 the situation changes, and the value σ 0 (µ) at which the potential V eff (σ; µ) is minimized acquires a non-trivial µ dependence.Therefore µ crit = σ 0 , the pole mass of the fermionic excitations in the zero-density theory.The value of the effective potential at the minimum acquires itself a µ dependence, so that the system supports a finite charge density.Let us start analyzing the small density limit.For µ > σ 0 and δµ ≡ µ − σ 0 ≪ σ 0 we find that the effective potential in a neighborhood of σ 0 is given by The potential is now minimized at so that the free energy density is given by Correspondingly, the charge density of fermions is given by

.26)
For small charge densities, we see that the discrete axial symmetry Z A 2 is in the broken phase.The symmetry breaking scale σ min (µ) decreases with increasing chemical potential and, correspondingly, charge density.This behavior can be understood by noticing that the finite µ contribution to V eff is a monotonic increasing function of σ, for positive σ, starting from a value of − N 12π 2 µ 4 for σ = 0 and increasing up to zero for σ ≥ µ (see figure 5).As a consequence, in the limit of large µ the Z A 2 symmetry will be restored, since the minimum of the effective potential V eff (σ; µ) will be in zero.This phase transition is in fact of first order and has an associated latent charge density, as can be readily verified by studying the full closed form effective potential we derived.A system with a charge density falling in between these two values is in a mixed phase and has chemical potential µ A .We provide a representative plot of the symmetry restoration phase transition in figure 6.
An analytic estimate for the critical value of chemical potential µ A corresponding to the Z A 2 symmetry restoration phase transition can be derived by comparing the value of the effective potential at the σ = 0 minimum, with the value of the µ = 0 effective potential at For our choice of parameters, the µ = 0 minimum is located at σ 0 ≈ 296, and we normalized the potential in such a way that its value at σ 0 is −1.In red: finite µ effective potential V eff (σ; µ), for increasing values of µ = 330, µ = 370 and µ = 400 (shown with growing intensity).There is a Z A 2 symmetry restoring first order phase transition at a critical value of µ, approximately given by µ A ≈ 370.For our choice of parameters, at the phase transition the charge density jumps from A , see eq. (5.30).
the minimum σ = σ 0 .Since the finite µ contribution is always non-positive, this estimate provides a lower bound on µ A .We find where in our estimate we have assumed σ 0 , and therefore N , to be large enough.The symmetry restoration chemical potential µ A is parametrically close to the critical chemical potential necessary to have a non-zero charge density, µ crit = σ 0 , with their ratio growing with σ 0 in a very mild way.For µ > µ A , the symmetry Z A 2 is restored, so that the fermionic excitations become effectively massless, and the properties of the system become equivalent to those of a system of N free massless Dirac fermions.Perhaps counterintuitively, the effect of the interactions induced by λ, g 4 , g 6 is overcome by the finite density contribution, driving the system to an effectively free finite density phase: (5.29) The charge density, in particular, is simply given by .30)We stress that this conclusion has general validity in this class of large N models, and holds even with the inclusion of higher order polynomial interactions for σ.Before concluding this section, let us comment on the question we started from: the possibility of U(1) V symmetry breaking at finite density.Our analysis of the Gross-Neveu model does not formally exclude this possibility, however it provides partial evidence for the hypothesis that the U(1) V symmetry remains unbroken.Indeed, the finite µ free energy has a non-trivial µ dependence and the system can support a finite charge density with unbroken U(1) V , differently from the scalar case analyzed in [1], where the breaking of U(1) V was shown to be necessary to support finite density in an O(N ) model.Moreover, the value of µ crit that we found coincides with the fermion pole mass, as expected on physical grounds, further hinting at the consistency of the scenario.Note that the argument discussed in sec.3.4 for Majorana order parameter can not be directly extended to the large N model, due to the U(N ) vectorial symmetry.Indeed, every Majorana-like fermion bilinear would transform in a non-trivial representation of the SU(N ) subgroup of U(N ), and cannot be regarded as an order parameter for U(1) V only.The simplest U(1) V scalar order parameter can be constructed, for N even, contracting the indices of N fermions with a completely antisymmetric tensor, and is never a fermion bilinear for N > 2.

Discussion and outlook
The dynamics of interacting systems of fermions at finite density and low temperature is far from being fully understood.A common belief is that, barring the case of free fermions (which form a Fermi gas), interacting fermions at low temperatures flow to a superfluid phase, or possibly to a non-Fermi liquid phase.The argument for this expected behavior relies on the existence of relevant deformations in the EFT of Fermi liquids [41][42][43], 6 that can generically lead to the formation of Cooper pairs, or more exotic strong coupling behavior.To the best of our knowledge, however, a complete understanding of the IR phases of interacting fermions is still missing, especially in the relativistic case, and the question of whether an interacting Fermi liquid can exist as a stable zero-temperature phase appears to be an open question.
In this work we described an approach to treat fermionic QFTs at finite chemical potential and zero temperature using path integral techniques.The method relies on an accurate treatment of the iε term needed to project on the ground state, and allows to compute finite µ quantities such as the free energy of a Fermi gas.We leveraged this technical tool to compute analytically the free energy of finite density QED in a homogeneous magnetic background, generalizing the Euler-Heisenberg effective action to finite density and reproducing the de Haas-van Alphen effect.These findings can be of interest in astrophysical contexts, such as strongly magnetized pulsars and relativistic plasmas.
As an application, we studied a generalization of the Gross-Neveu model to 3 + 1 spacetime dimensions, and showed that in the large N limit it can support a finite density phase with unbroken U(1) V internal symmetry.This is to be contrasted with the case of the scalar O(N ) vector model, where a finite density phase can only exist along with the spontaneous breaking of the U(1) charge associated to µ [1], leading to a superfluid behavior.This suggests that there might exist theories of interacting fermions, in physical spacetime dimensions, with a zero-temperature Fermi liquid phase.A definite answer requires a stability analysis, together with a study of 1/N corrections and the dynamics of excitations.We hope to address these questions in the future, and make further progress towards the goal of a general understanding of the infrared phases of relativistic quantum fields at finite density.where q 2 max = 2eB(ℓ max − 1), and we used the fact the F (ℓ max ) ∼ O (eB) 5/2 and F (p) (ℓ max − 1) ∼ O (eB) 5/2 for p = 0, 1, 2. One can now recognize the integral as an integral in cylindrical coordinates over a domain Ω which is the intersection of a ball B µ of radius µ 2 − m 2 and a concentric cylinder C B of infinite length and radius q max .The boundary of Ω is the Fermi surface ∂Ω.We can write the free energy as (B.3) In the limit 2eB → 0 the radius of the cylinder tends to q max = µ 2 − m 2 and Ω tends to the ball B µ .We thus recover the B = 0 result of equation (3.24), up to a µ-independent cosmological constant term coming from f (B = 0, 0) that we drop.For small non-zero magnetic field we approximate the sum as an integral and write the result as a correction to the B = 0 case f (B, µ) max is a small parameter in the limit 2eB → 0, of order ∆ ∼ O( √ eB).We can expand this result in powers of ∆ obtaining

4 10 45
Absence of symmetry breaking for free fermions The Fermi gas in a magnetic background 12 4.1 The functional determinant 12 4.2 A modern derivation of the nonperturbative Euler-Heisenberg effective action 13 4.3 Finite density QED in a magnetic background 16 4.4 Magnetic susceptibility and the de Haas-van Alphen effect 18 The four dimensional Gross-Neveu model at finite density 20 5

Figure 2 :
Figure 2: Magnetic susceptibility of the zero density QED vacuum in an external magnetic field, in units of α/π and as a function of β ≡ 2eB/m 2 .The dashed red (N = 1) and green (N = 2) lines provide truncations of the small β asymptotic expansion with N terms.Increasing N improves the accuracy for small β, but reduces the region in which the approximation is valid.The strong field approximation (shown as an orange dot-dashed line) provides an accurate result only for β ≳ 100.The intermediate region is described by the full non-perturbative result.

Figure 3 :
Figure 3: Magnetic susceptibility of the finite density QED vacuum in a strong external magnetic field, for different values of μ = µ/m.We use variables that make manifest the filling of the first five Landau levels.The magnetic susceptibility has a periodic spike feature as a function of 1/β, occurring for integer filling, known as the de Haas-van Alphen effect.Continuous lines: analytic form of the finite µ contribution (second term in eq.(4.47)).Dots: full result, evaluated at equidistant points.

Figure 4 :
Figure 4: Same as figure3but for a larger range of 1/β and for μ2 = µ 2 /m 2 = 3/2.In the weak field limit (small β), the magnetic susceptibility χ B oscillates more and more wildly.The smooth result (4.51) is obtained only as an average.

Figure 5 :
Figure 5: Finite µ contribution to the effective potential of σ in the (generalized) Gross-Neveu model in 3 + 1 dimensions, eq.(5.20).The result is normalized in units of N µ 4 and expressed as a function of σ/µ.

Figure 6 :
Figure6: In blue: effective potential for a (generalized) Gross-Neveu model, in units M = 1, for κ 2 = 1, κ 4 = 2, κ 6 = 0.5 and N = 10 6 , see eq.(5.13).For our choice of parameters, the µ = 0 minimum is located at σ 0 ≈ 296, and we normalized the potential in such a way that its value at σ 0 is −1.In red: finite µ effective potential V eff (σ; µ), for increasing values of µ = 330, µ = 370 and µ = 400 (shown with growing intensity).There is a Z A 2 symmetry restoring first order phase transition at a critical value of µ, approximately given by µ A ≈ 370.For our choice of parameters, at the phase transition the charge density jumps from Q − (µ A ) ≈ 0.25 to Q + (µ A ) = 1 in units of