Beyond N =∞ in Large N Conformal Vector Models at Finite Temperature

We investigate finite-temperature observables in three-dimensional large N critical vector models taking into account the effects suppressed by 1 N . Such subleading contributions are captured by the fluctuations of the Hubbard-Stratonovich auxiliary field which need to be handled with care due to a subtle divergence structure which we clarify. The examples we consider include the scalar O(N) model, the Gross-Neveu model, the Nambu-Jona-Lasinio model and the massless Chern-Simons Quantum Electrodynamics. We present explicit results for the free energy density to the subleading order, which also captures the one-point function of the stress-energy tensor, and include the dependence on a chemical potential. We further provide a formula from diagrammatics for the one-point functions of general singletrace higher-spin currents. We observe that in most cases considered, these subleading effects lift the apparent degeneracies between observables in different models at infinite N , while in special cases the discrepancies only start to appear at the next-to-subleading order. September 6, 2023 ar X iv :2 30 9. 02 34 7v 1 [ he pth ] 5 S ep 2 02 3


Introduction and Summary
The understanding of strongly coupled quantum matter at finite temperature is an important problem in physics. Its broad applications range from the basic experimental needs, since any realistic quantum critical point would necessarily be at non-zero temperature, to fundamental theoretical questions such as the thermalization of many-body quantum systems, where highly excited states of the system are conjectured to be universally approximated by the finite-temperature (thermal) state [1,2], and furthermore to profound connections between quantum gravity and quantum field theory dictated by the holographic principle [3][4][5], where finite-temperature systems are believed to be dual to black holes, and therefore would provide valuable information about the quantum nature of the latter.
Conformal Field Theory (CFT), which provides a non-perturbative formulation of critical phenomena, offers a powerful approach to investigate the finite-temperature quantum system near its quantum critical point. While the second-order quantum phase transition occurs at zero temperature and is described by the CFT on the flat spacetime R d−1,1 , to turn on finite temperature T amounts to compactifying the Euclidean time τ ∼ τ + β with β = 1 T and imposing suitable periodicity conditions along τ . At thermal equilibrium, this boils down to studying the Euclidean CFT (from Wick rotation) on S 1 β × R d−1 (known as the thermal background), and the basic observables are correlation functions of local operators on this background, the simplest of which being the one-point function ⟨T µν ⟩ β of the stress-energy tensor which measures the free energy of the thermal state. On the one hand, such thermal observables behave qualitatively different from the conventional flat space CFT correlation functions due to the explicit breaking of conformal symmetry by the thermal background, and for the same reason they appear much harder to determine in interacting models. On the other hand, these thermal observables present a universal coarse-grained description of the flat space operator data, through the asymptotic density of states in the CFT Hilbert space on S d−1 and operator-product-expansion (OPE) coefficients averaged over high energy states [6] (see also the recent work [7]). In light of the AdS/CFT correspondence [3][4][5], this translates to the expectation that the black hole solution in gravity is a universal coarsegrained description that encodes a large number of underlying microstates. Nonetheless deriving such universal formulas directly from the flat space OPE data is very challenging in interacting CFTs of spacetime dimension d ≥ 3. A main focus of this work is to provide explicit solutions to thermal observables using field theory techniques in interacting CFTs in d = 3. 1 1 It would be interesting to study our results in relation to the flat space OPE data. We leave that to There is an alternative interpretation of the CFT on S 1 β × R d−1 entirely as a zerotemperature quantum system (or a classical statistical system in d dimensions) and consequently a reinterpretation of the results we present in this paper. Instead of taking S 1 β to be the compactified Euclidean time, we can choose it to be a compactified spatial circle (and keep the new Euclidean time along one of the R d−1 directions). This way, it describes the quantum system confined to a finite interval of length β with certain periodic boundary conditions, also known as the Kaluza-Klein compactification of the CFT. Here an important and universal observable is the Casimir force, which underpins the quantum nature of the system. For instance, if we confine free electromagnetic field in cavity, the walls of the cavity would feel the pressure due to the quantum fluctuations of the vacuum despite being at zero temperature. Such an effect was firstly predicted by Hendrik Casimir in 1948 and later confirmed experimentally. The Casimir effect exists for any quantum field theory, but in the absence of a mass-gap, the Casimir pressure is expected to depend on the geometric moduli of the cavity algebraically. Consequently the Casimir force will be long-ranged and more easily detected in a gapless phase, in contrast to a gapped phase where the Casimir pressure decays exponentially as we increase the cavity size. Therefore, the Casimir force (pressure) serves as a salient observable when the system undergoes a quantum phase transition (or a second-order phase transition in the statistical model), which is often described by a CFT.
It is again measured by the one-point function ⟨T µν ⟩ β but with the time and space directions swapped. Similar reinterpretations hold for more general observables on S 1 β × R d−1 . In this work we study d = 3 CFTs on the background S 1 β × R 2 . The dimensionful parameter β breaks the conformal symmetry explicitly to the Euclidean isometry on R 2 and translation symmetry along S 1 β . As an immediate consequence, the one-point functions of local operators are no longer constrained to be zero. Instead the residual symmetry and dimensional analysis require the one-point function of a primary operator O µ 1 ...µ ℓ of scaling dimension ∆ and spin ℓ to take the following form with an overall constant b O [6], For instance, the stress-energy tensor has the following one-point function at finite temfuture work.
perature (we keep d general for the moment), where the constant f = b T d determines the free energy density of the CFT thermal state via which is negative by the positivity of energy [6]. 2 When the S 1 β is regarded as a spatial circle, f is referred as the critical Casimir amplitude which determines the Casimir pressure from changing the circle size β. This universal quantity f have been measured for a variety of systems in experiments or from the Monte-Carlo simulations (see [8] for an extensive review).
Yet only limited results are available from the theoretical side for d ≥ 3 [6,[9][10][11][12][13][14][15]. One main purpose of this paper is to compute the free energy coefficient (equivalently critical Casimir amplitude) f explicitly in interacting d = 3 CFTs that are solvable in some regime, and similarly for one-point functions (1.1) of more general local operators.
Importantly these one-point functions also give access to flat space CFT data that are hard to obtain otherwise. Regarding the background S 1 β × R d−1 as a limit of S 1 β × S d−1 R as R β → ∞, we can see the free energy coefficient controls the asymptotic density of high energy CFT states on S d−1 , equivalently heavy local operators by the state-operator correspondence.
Explicitly, the partition function S 1 β × S d−1 R that counts states in the Hilbert space H S d−1 graded by dilatation operator ∆ is determined by the free energy coefficient f in this limit, where ρ(∆) is the density of states and V d−1 = 2π R d−1 is the area of S d−1 R . Performing an inverse Laplace transform, this implies the following asymptotic density of heavy operators [16], by the thermal one-point function ⟨O⟩ β for O ≡ e µ 1 . . . e µ J O µ 1 ...µ J in the same limit, where ϕ labels an orthonormal basis in H S d−1 . The inverse Laplace transform then produces an asymptotic formula for the averaged OPE coefficients [17], Note that all of the above are completely determined by the one-point function coefficient b O (including the free energy coefficient f as a special case).
In general, CFTs are strongly coupled and thus inaccessible via small perturbations near some integrable theory. A class of d = 3 models with a large number N of scalars and fermions governed by O(N ) (or U (N )) invariant interactions, which we will refer to collectively as large N vector models, circumvent this obstacle by admitting an expansion in 1 N and thus providing an ideal playground to study interacting CFTs (see [18] for a review). Surprisingly, this simple-looking vector model has proven to be an excellent description of real-world critical systems. In the case of scalar vector models, the UV Lagrangian for the scalar fields ϕ i with i = 1, 2, . . . , N takes the following form The special cases with N = 1, 2, 3 correspond to the Ising model, the XY model and the Heisenberg model respectively, each of which has experimental realizations both as zerotemperature quantum phase transitions and classical (thermal) phases transitions in statistical systems and can also be simulated on a lattice. 3 Furthermore, the limit N → 0 describes the statistics of polymers [20]. In the case of the Ising and the XY model, the free energy 3) (equivalently the critical Casimir amplitude) has been computed by the Monte-Carlo techniques [21,22] and using the ϵ expansion and functional renormalization group [23]. The results are summarized below, In the large N limit, the leading contribution to the coefficient f was determined in [24], Famously, the 1 N expansion in the large N scalar vector model (similarly for fermionic models) coincides with a semi-classical expansion where 1 N plays the role of the Planck constant ℏ [18]. This semi-classical expansion is facilitated by introducing the Hubbard-Stratonovich auxiliary field σ and rewriting the action (1.8) in a form that is quadratic in the elementary fields ϕ i . Path-integrating over ϕ i then produces a non-local effective Lagrangian F(σ) for the σ field which is amenable to saddle-point analysis due to the small ℏ ∼ 1 N . In the scaling region, the saddle-point σ = σ * in the N → ∞ limit and the value of the effective Lagrangian F(σ * ) determines the O(N ) CFT completely to the leading order in N . The first correction comes from the one-loop determinant around the saddle σ = σ * , which is suppressed by 1 N . For the large N critical O(N ) scalar model, the 1 N correction to the free energy coefficient f in (1.3) was first computed in [9]. The calculation requires a subtle numerical procedure which we will clarify here. We also provide several approaches to access the same observable.
Together they lead to the following expression for the free energy coefficient for the O(N ) model including the first 1 N correction, which improves significantly on the precision of the result from [9]. It is immediate to note that the correction is comparable to the leading term at small N . This is related to the fact that we solve these models using the saddle-point approximation, which generally produces asymptotic series in the expansion parameter.
We further analyze the 1 N correction to the free energies of fermionic large N vector models using our numerical procedure, including the Gross-Neveu (GN) model and the Nambu-Jona-Lasinio (NJL) model of N Dirac fermions with different four-fermion interactions. The results are summarized below, where Cl 2 is the second Clausen function (see (3.17)). We see that the two models has identical free energy density at leading order in the large N expansion whereas the 1 N correction tells them apart. We have also analyzed the same observable for the chiral Ising GN model, which is a close cousin of the GN model with another different four-fermion coupling and reduced global symmetry. Yet its free energy density only starts to differ from the GN model at the next-to-subleading order in 1 N . The fermionic large N models present additional subtleties that are absent in the bosonic models. For instance, the large N saddle-point equation in the scaling region has additional spurious solutions that do not describe the actual CFT. We observe that this can be diagonosed by studying the 1 N corrections and only in the correct saddle-point (corresponding to the CFT) the dependence on the UV regularization vanishes, as expected for a scale invariant theory.
We also study an interesting class of 3d CFTs that share the same U (N ) global symmetry.
They are described by the 3d Chern-Simons Quantum Electrodynamics (CSQED) with N charge one Dirac fermions and Chern-Simons level k. In the 't Hooft limit, namely N, |k| → ∞ with the 't Hooft coupling λ = 4πN k fixed, they give rise to a one-parameter family of CFTs labelled by λ. We find that while the free energy coefficient in this case is independent of λ in the leading large N order, the subleading piece is a nontrivial function g(λ) as in, which has the following limiting behaviors, 15) corresponding to QED with vanishing Chern-Simons level k = 0 and the free limit (with infinity Chern-Simons level) respectively. The former agrees with the previous result in [12].
In addition to the free energy coefficient f in (1.3), we also develop diagrammatic methods in the large N vector models to compute the one-point functions of other primary operators as in (1.1). Because of the unbroken global symmetry, only singlet operators with respect to the symmetry acquire nontrivial one-point functions. In the large N vector models, a family of such operators are known as the single-trace higher-spin currents J s µ 1 ...µs which takes the following schematic form in the scalar vector model [25][26][27][28][29][30], with even spin s. These operators are conserved currents in the N = ∞ limit but develop anomalous dimensions at the order 1 N . Their thermal one-point function coefficients b s in the N = ∞ limit of the scalar O(N ) model has been previously computed in [6] using the inversion formula. Here we provide a direct diagrammatic derivation of b s for both bosonic and fermionic large N vector models in the leading large N limit. We also describe explicitly the procedure to obtain their 1 N corrections where the inversion formula in [6] does not obviously apply.
We now discuss a number of potential applications and future directions of our work.
As the Monte-Carlo simulations are being pushed to study vector models at higher N , it would be interesting to see how our results on the 1 N corrections compare to these numerical investigations. It would also be interesting to see if these large N systems can be realized in experiments (e.g. by stacking and twisting existing finite N setups), which would measure the 1 N corrections we have found. Another interesting direction involves treating the parameter N as a coupling constant and finding the non-perturbative corrections to the large N results using the techniques developed by Lipatov [31,32] and renormalons [33]. As special instances of the AdS/CFT correspondence [3][4][5], the d = 3 large N vector models are expected to be dual to certain versions of Vasiliev's higher-spin gravity on AdS 4 [34][35][36][37][38] (see also [39,40] for recent works). The critical large N vector model at finite temperature should then correspond to a black brane solution of Vasiliev's higher-spin gravity. Such a solution is known in the linearized limit, but not yet in the full non-linear theory due to the high complexity of the bulk interactions in the higher-spin gauge theory [41,42]. Furthermore, it is not known how to reproduce even the leading N free energy from a bulk action for the higher-spin gravity. 4 Nonetheless in light of the successful matching between one-loop effects (i.e. order O(N 0 ) contributions to the free energy) for the CFT on S 3 and the higher-spin gravity on AdS 4 in [49], it would be interesting to compare the one-loop contributions from higher-spin gauge fields on a putative thermal geometry and the finite temperature 1 N corrections we find on the CFT side, which will provide a nontrivial test on the bulk solution and offer some insights on the structure of the bulk solution. Finally, it was shown recently in [50] that by incorporating global symmetry twists one can obtain a refined version of the asymptotic density of states (1.5) labelled the representations of the symmetry. Here we consider the simplest twist by turning on the chemical potential for a particular U (1) symmetry of the CFT. It would be interesting to study such symmetry refinement in a more general setting and also compare with explicit operator counting.
The rest of the paper is organized in the following way. In Section 2, we consider the critical scalar O(N ) model at finite temperature in detail. We compute the two-point function of the σ field at finite temperature, and then numerically obtain the first correction to the free energy coefficient in the large N expansion and comment on subtleties in the numerical procedure. We also include the dependence of the free energy on a particular U (1) chemical potential. In addition, we show explicitly how the same result follows from the computation of the one-point function of the stress-energy tensor by explicitly summing Feynman diagrams. We then generalize the diagrammatic analysis to higher-spin currents of even spin s > 2. In Section 3, we apply the developed techniques to study fermionic vector models including the Gross-Neveu model and the Nambu-Lasino-Jonas model. In Section 4, where m 0 and λ 0 are the bare mass and coupling and ϕ i belongs to the vector representation of the O(N ) global symmetry group. Famously, the model admits a large N limit where physical observables such as the correlation functions of local operators can be extracted using a saddle-point approximation (see [18] for an extensive review). To see that, we apply the Hubbard-Stratonovich (HS) trick using the auxiliary σ field, also known as the HS field, whose original integration contour is along the imaginary axis [18]. Then up to a constant shift, we have with r 0 = m 2 0 λ 0 . Note that σ field can be thought as a mass of the field ϕ i . Integrating out the fields ϕ i , we arrive at the following effective Lagrangian for σ, By demanding λ 0 = λ t 0 N , r 0 = N r t 0 with λ t 0 , r t 0 held fixed in the large N limit, we see that 1 N plays the role of an effective Planck constant which gives rise to a semi-classical expansion.
For that purpose we first need to solve the equation of motion for σ, where the LHS coincide with the coincident limit of the propagator for ϕ i , Assuming the homogeneous ansatz σ(x) = σ, the saddle-point equation for σ is given by, which is also known as the gap equation since σ determines the mass gap for the scalar fields. The LHS of the above equation is divergent for d ≥ 2 and needs to be regularized.
It is straightforward to check that the divergences could be absorbed in the redefinition of the bare coupling constants r t 0 . To study the CFT, we bring the system to a critical point by further fine-tuning this parameter. Note that for d < 4 near the free Gaussian point, the composite operators ϕ i ϕ i and (ϕ i ϕ i ) 2 are both relevant and thus we should fine-tune them at the same time to reach the free fixed point. Perturbed away from the Gaussian point by the (ϕ i ϕ i ) 2 operator, the theory can flow to an interacting critical point, where (ϕ i ϕ i ) 2 becomes irrelevant. In this case, we only need to fine-tune the mass r 0 . This critical point describes the second-order phase transition between the ordered and the disordered phases of the O(N ) model. The actual value of the parameter r 0 where the system becomes critical is scheme-dependent and we will label the regularization scheme by R. The phase transition occurs at σ = 0 (for σ < 0 there are tachyonic instabilities) since σ controls the correlation length. It is thus more convenient to parameterize the bare coupling r t 0 as where the RHS is computed with the chosen regularization scheme in the UV. Note that this equation has a solution only for d > 2 since at d = 2 the integral acquires an IR divergence.
This is a manifestation of the Coleman-Mermin-Wagner theorem stating that there are no Goldstone modes in d ≤ 2. Substituting the value of r t 0 (2.7) in the gap equation (2.6), we obtainˆR Now by tuning r t we can bring the system to the critical phase transition. Indeed, with σ > 0, we obtain the following equation .
Expanding the RHS of the above equation for small σ at 2 < d < 4 gives [18] where K d is a scheme-independent constant and I R is a constant that depends on the regularization scheme R. The terms that are subleading when σ → 0 are omitted in (2.10). Near the critical point, we have the scaling If we further set λ t 0 = − 1 I R then the corrections due to the finite-size effects would become suppressed. This value should be carefully chosen when adopting a numerical lattice regularization scheme (see Appendix A for further comments). For the analytical treatment, we find the dimensional regularization to be the most convenient. For so to reach the critical point we set r t = 0 (equivalently r t 0 = 0) and to cancel subleading corrections we pick λ t 0 = ∞ in this scheme. Consequently, the corresponding CFT on flat space is governed by the following action and similarly the effective Lagrangian for σ at the critical point is (2.14) To determine observables at the subleading order in the 1 N expansion, we also need the propagator for σ, which follows from the large N Lagrangian (2.3), The negative signs reflect the fact that the original integration contour for σ runs parallel to the imaginary axis [18]. At the critical point in the dimensional regularization scheme (see (2.12)), we have from (2.5), so that the scaling dimension of the composite operator σ ∝ (ϕ i ) 2 is ∆ σ = 2 and independent of the spacetime dimension (thus different from the mean-field value for general d, a hallmark of interacting CFT).
Here we are mostly interested in the d = 3 CFT on the background S 1 β × R 2 which is flat but has nontrivial topology. This is described by the same action (2.13) with periodic boundary conditions for the fields along S 1 β . The gap equation (saddle-point equation) for σ now takes the following form, 5 where ω n ≡ 2πn β are the bosonic Matsubara frequencies and the UV divergence is regulated in the same way as in flat space. As we review in the next subsection, the saddle-point solution σ = σ * is no longer zero for finite β and corresponds to a finite mass-gap generated by quantum effects. Equivalently, it determines the one-point function of σ ∼ (ϕ i ) 2 in the large N limit, where (ϕ i ) 2 ren is the renormalized primary operator with the normalized two-point function on flat space (2.19)

Large N Free Energy and the Subleading Correction
Here we derive the free energy density of the scalar O(N ) CFT on S 1 β × R 2 , where V 2 regulates the infinite spatial volume and we focus on the non-extensive (in S 1 β ) part of log Z S 1 β ×R 2 which is free from counterterm ambiguities. The free energy density has the following large N expansion, where F 0 denotes the leading contribution of order O(N ) and F −1 is the first subleading contribution at O(1).
We start with the effective Lagrangian for σ,

22)
5 In the following we will simply write ⃗ k as k for the two dimensional momentum to simplify the notation.
where ω n ≡ 2πn β . In the large N limit, the leading free energy density is determined by the saddle point σ = σ * , where we work in a general regularization scheme R, keep only terms that are non-vanishing for constant σ (sufficient for the leading large N analysis) and restore r t 0 (which vanishes in the dimensional regularization). The last term in the above expression is the cosmological constant counterterm.
A key point is that the sum-integral in (2.24) has the following large energy momentum expansion, where the omitted terms are absolutely convergent thus independent of the regularization scheme R. Furthermore, the first term in the bracket is also absolutely convergent which follows from (2.6) at σ = 0 (zero temperature). Consequently F ren (σ) is scheme-independent.
To facilitate the analytic calculations, we perform dimensional regularization on the twodimensional momentum integral and zeta function regularization on the Matsubara sum.
Implementing this procedure to (2.22), we first obtain The following identity (which holds in zeta function regularization) will be useful [51], Integrating in σ, we then obtain the renormalized Lagrangian 6 The integration constant is fixed by requiring, which follows directly from (2.26) by zeta function regularization.
Extremizing F ren (σ) with respect to σ, we find the thermal gap equation (the regularized version of (2.17)) with the following solution Consequently, the leading free energy density for the O(N ) CFT is Let us briefly comment on a subtle point in this calculation. The CFT is Lorentz invariant at zero temperature (i.e. β = ∞). It is a priori not obvious whether the short-distance regulator implemented above (a combination of zeta function and dimensional regularizations) respects the Lorentz symmetry at scales much smaller than β. In fact, as we will see momentarily, this regulator is not compatible with the Lorentz symmetry at the subleading order in the 1 N expansion. Nonetheless, as explained after (2.24), in the leading large N limit, the scheme independence is enhanced.
As a consistency check, we also evaluate (2.24) with a manifestly Lorentz-invariant regulator [9] by performing sum-integral with the following Lorentz-invariant hard cutoff Note that in a general regularization scheme we will need to subtract off the extensive piece of the free energy to obtain the scheme-independent part of the free energy density F (see (2.20)). Such an extensive contribution to F in a CFT on S 1 β × R 2 can only come from the cosmological constant (β independent). For our choice of regularization scheme here, such a constant is absent, since the last two terms in (2.24) vanish in this scheme. and obtain the same results numerically as in (2.31). 7 We are mostly interested in the first subleading correction in 1 N denoted by F −1 (β) in (2.21). According to the semi-classical expansion in the large N vector model, this is computed by the log-determinant of the second variation of the effective Lagrangian (2.3) with respect to σ (equivalently from the inverse propagator G −1 σ ), subject to regularization and renormalization that we explain below, and Π β takes the following form with σ * as in (2.31). We take the integral over the spatial momentum using the Feynman parametrization and evaluate the sum over the discrete frequencies exactly. The detailed computation is presented in Appendix B, and the resulting expression is To evaluate (2.34) using (2.37) requires a further regularization. Indeed, we see that the σ self-energy at large momentum behaves as (see Appendix B.1 for details), where The first term in this expansion (2.38) corresponds to the flat space time propagator and 7 We emphasize that r t 0 ̸ = 0 at the critical point in this regularization scheme.
should be cancelled to obtained the renormalized free energy. We thus arrive at the following expression for the renormalized free energy at the subleading order in 1 N , where the last term in the second line follows from (2.29) (see Appendix B for further details).
The second term in (2.38) is proportional to the gap equation (2.30) and vanishes for the special point σ = σ * from (2.31) relevant for describing the finite temperature CFT. If this term were nonzero, we would encounter a quadratic divergence in the subleading free energy (2.40) that is linear in the temperature, which would be inconsistent with the structure of UV divergences in local quantum field theory.
The next term in the large momentum expansion (2.38) no longer vanishes at the critical point σ = σ * , and consequently The further subleading terms in the above expansion converge absolutely and we don't need to worry about them. On the other hand, the first term on the RHS is dangerous since it contributes an apparent logarithmic divergence to (2.40) and could lead to regularization ambiguities. To see this explicitly, we consider the following sum-integral, procedure and consequently should use a Lorentz-invariant UV regulator, such as (2.33), which corresponds to choosing M n = {p|p 2 + ω 2 n < Λ 2 } in (2.42). 8 This gives which determines the contribution from this term to the free energy (2.40). In the next section, we will evaluate (2.40) numerically taking into account this regulator subtlety.

Numerical Calculation
Now we are in position to compute numerically the subleading correction to the free energy for large N critical vector model, by evaluating the sum-integral in (2.40). Because of the oscillatory behavior of Π β (see (2.37)), this sum-integral needs to be handled with care. Below we explain the strategy for this numerical evaluation which we implement in mathematica.
We start by rescaling all momentum and energy in (2.40) by β, such that the β dependence is completely factored out and given by 1 β 3 . The goal is to determine the dimensionless coefficient in F −1 (β).
We first need an efficient way to evaluate the σ self-energy Π β (ω n , p) to high precision.
The expression (2.37) for Π β (ω n , p) is a highly oscillating integral at large momentum, thus converges very slowly. Therefore, to optimize the numerical evaluation of Π β (ω n , p), we introduce two spherical shells of radius Λ 0 and Λ, with Λ 0 < Λ. Where Λ is the UV cut-off and Λ 0 is chosen in such a way that the difference between the actual value of the self-energy (2.37) and its large momentum expansion (B.10) would be negligible.
Between the two spherical shells Λ 2 > ω 2 n + p 2 > Λ 2 0 , we can use the large P expansion of the σ self-energy to compute the sum-integral to the desired precision, by keeping all terms in (B.10). Finally, the relativistic UV cutoff Λ is taken to infinity numerically. Note that this resolves the regularization ambiguity which we have discussed near the end of the last section.
Implementing the procedure outlined above in Mathematica and taking different values of 250 ≤ Λ 0 ≤ 450 and 10 4 ≤ Λ ≤ 10 8 , we find that the sum-integral in (2.40) evaluates to and consequently the subleading correction to the free energy of the O(N ) vector model reads,

Turning on Chemical Potential
One natural extension of our analysis of the finite temperature free energy of the O(N ) CFT in the previous section is to include a background for its global symmetry. Here for simplicity, we consider N even and take the U (1) subgroup of O(N ) with commutant SU (N/2) (such that the complex scalars ϕ j + iϕ j+N/2 for j = 1, . . . , N/2 have charge 1 under this U (1)). We turn on an imaginary chemical potential parameterized by µ ∈ [0, 2π) for this U (1) subgroup, via a background gauge field with nonzero temporal component A 0 = µi β . Consequently the energy spectrum for the scalar fields are shifted to βω n = 2πn + µ. The effective Lagrangian in this case is and its renormalized version reads (2.48) When µ = 0, this reduces to the case considered previously (see (2.28)). The gap equation is given by, We denote its solution byσ =σ * (µ) whose explicit form is given below, The free energy in the leading large N limit follows from which is plotted in Figure 1. 9 Moving onto the subleading order in the 1 N expansion, we will determine the free energy F −1 (β, µ) by performing the sum-integral as in (2.24) which now involves the σ free-energy Π µ β that depends on µ (see (B.4) for its explicit integral representation). Implementing the numerical procedure explained in Section 2.3, we compute the F −1 (β, µ) as a function of the chemical potential and the result is presented in Figure 2.

One-point Function of the Stress-energy Tensor
As explained in the introduction, the thermal free energy of the CFT contains the same information as the one-point function of the stress-energy tensor at finite temperature (see around (1.2)). So far we have analyzed the free energy of the O(N ) CFT using the semi-classical expansion in the large N limit. It would be useful to understand how the same physical quantity can be computed using the standard Feynman diagrams for this Lagrangian field theory. In particular, this second diagramatic approach will have immediate generalizations to determining the thermal one-point functions (1.1) of more general operators in the CFT, which we will discuss in the next section.
T 00 T 00 Figure 3: Examples of Feynman diagrams (on the right) containing ϕ i loops that contribute to the thermal one-point function ⟨T 00 ⟩ 0,β at the leading order in 1 N . They resum to the diagram on the left with the large N exact propagator for ϕ i represented by a thick line.
T 00 T 00 T 00 T 00 Figure 4: Examples for the four families of Feynman diagrams that can contribute to the thermal one-point function ⟨T 00 ⟩ −1,β at the first subleading order in 1 N . All internal lines are large N exact propagators for ϕ i .
From the general structure of the thermal one-point function (1.2), it suffices to focus on the temporal component T 00 of the stress tensor T µν , which in the O(N ) CFT is given by, up to the improvement total derivative terms, 10 The explicit relation between the thermal one-point function and the free energy in the d = 3 CFT is and similar to the RHS analyzed previously (2.21), the LHS admits an 1 which can be seen explicitly by reorganizing the Feynman diagrams at large N .
The leading large N result ⟨T 00 ⟩ 0,β for the thermal one-point function only receives contributions from the expectation values of the first two terms in (2.52) and comes from resummation of the Feynman diagrams with increasing number of ϕ i loops (see Figure 3), which induces a non-zero value for the sigma field. Therefore we have, from the one-loop diagram with the large N exact ϕ i propagator, which confirms the relation (2.53) upon comparing with (2.22) with σ = σ * as in (2.31).
We now compute the 1 N correction to the thermal one-point function ⟨T 00 ⟩ β . The Feynman diagrams that contribute at this order come in four families which are given in Figure 4 with increasing number of ϕ i loops with exact propagators, independent of the temperature (or more general spacetime background). Each family of these diagrams can be resummed in the IR using the σ self-energy (inverse propagator) Π β (ω n , p) (which captures the contributions from a chain of ϕ i loops) given explicitly in (2.37) and together they determine ⟨T 00 ⟩ −1,β , (2.56) Explicitly the first two contributions could be written as The second equation above can be further simplified to where we have used in the last equality that the saddle σ = σ * in (2.31) satisfies The contributions from the first two terms in (2.56) then combine to, (2.60) Similarly, the third family of diagrams in Figure 4 resum to Finally, the fourth family of diagrams in Figure 4 vanish, where G reg (x, x) is the regularized two-point function of ϕ i at coinciding points which vanishes as a consequence of the gap equation (2.17). Combing the above all together, we arrive at the full subleading 1 N correction to the stress tensor one-point function, (2.63) Let us compare the above expression with that of the free energy. We have from Section 2.2 where the last term clearly coincides with (2.63) after using (2.36).

One-point Function of Higher-Spin Currents
Using the ideas from the previous section, we can compute the 1 N corrections to the thermal one-point function of more general local operators, which, as emphasized in the introduc- up to total derivatives which are fixed by demanding the LHS to be a primary operator.
They generalize the stress-energy tensor which appears at s = 2 (as well as the singlet scalar operator (ϕ i ) 2 ren at s = 0). Using the residual symmetries of the thermal background, the one-point function of the higher-spin current of rank s is constrained to take the following form where γ s is the anomalous dimension which is suppressed by 1 N in the large N limit and thus the O(N ) CFT is said to have a slightly broken higher-spin symmetry [27]. It is convenient to introduce an auxiliary complex null polarization vector ξ µ = ξ (1, i, 0) and write 12 whose one-point function is where we have used (2.67). Therefore to determine the one-point function of the higher-spin current J s µ 1 ...µs , it suffices to focus on the component J s ξ . In the large N limit, the thermal one-point function ⟨J s ξ ⟩ β admits a 1 N expansion similar to (2.54) for the stress tensor. At the leading order, the same family of diagrams as in Figure 3 contribute and give where σ * is given in (2.31). As usual, the UV divergences in the sum-integral above is regulated by subtracting the flat space expression in the given regularization scheme. To 12 Note that these are bare operators (also (2.66)) which will be subjected to renormalizations (suppressed by 1 N  simplify the computation, we introduce the generating function where the regularization of the last term above is implicit. This generating function is just a specialization of the coordinate-space propagator for ϕ i , which for general any spacetime separation r µ = (r 0 , ⃗ r) is given by which follows from (2.72) by Poisson resummation. Physically, each summand above represents the contribution from a worldline instanton of mass √ σ * = ∆ β propagating along a geodesic on the thermal background that connects the two ϕ i insertions. The m = 0 term coincides with the flat space-time propagator in the short distance limit. After subtracting that, we obtain the regularized version of G β,0 (ξ), which gives The one-point function coefficient b s in the leading large N limit (which we will denote as b 0 s ) of the spin s current is obtained immediately by the small ξ expansion of the RHS above, which agrees with the results obtained in [6] from a different method (up to an overall normalization factor that only depends on s).
Alternatively, we can directly regulate the integrals (2.71) as follows (setting β = 1 below and using (2.31)), By standard contour manipulation, the above can be rewritten as where ϵ p ≡ p 2 + ∆ 2 and γ is a contour that circles the simple poles at z = 2πin with n ∈ Z in the counter-clockwise direction. 13 Deforming the contour to infinity, we obtain by for positive even s, which again reproduces the same result (2.75), after implementing a straightforward change of variables and using the following integral identity The first subleading correction to the one-point function of higher-spin currents The subleading correction to the one-point function of higher-spin currents can be computed using Feynman diagrams as in the case of the stress-energy tensor in Section 2.5.
The dynamical data in the one-point function (2.70) of the higher-spin currents, namely the overall constant coefficient b s and the anomalous dimension γ s (see (2.68)) have the The 1 N corrections ∆ 1 s to the scaling dimensions of these operators have been computed in [25,29,30,52], Consequently, the first correction in the 1 N expansion of the one-point function of bare operator (2.70) takes the following form, where Λ is a UV regulator and the UV divergence could be treated by proper renormalization of the operator. In the above equation, we have reintroduced the subscript on the higher-spin current J s ξ,bare to emphasize this is the bare operator. For the stress-energy tensor, its scaling dimension does not receive corrections (∆ 1 2 = 0) and therefore the logarithmic correction is absent, which simplifies the computation. 14 Following the strategy of the previous section, we start by incorporating the first 1 N correction to the ϕ i propagator as follows, where Σ(ω n , p) is a self-energy part for the matter field ϕ i where σ * is the large N saddle point for σ (as given in (2.31) in the dimensional regularization). Note that the self-energy Σ(ω n , p) also contains the information about 1 N corrections to the one-point function of σ ∼ (ϕ i ) 2 (i.e. correction to σ * ).
The bare self-energy Σ(ω n , p) is computed from the second and third diagrams in the Figure 5, which produce Then the one-point function to this order in the 1 N expansion is then Note that the above expression contains a UV divergence that renormalizes the operator itself, which could be accounted for by the following counter-term to the first subleading where ∆ 1 s is given in (2.81). This incorporates contributions from all diagrams in the Figure 5. We then obtain the following sum-integral expression for the correction to the one-point .

(2.88)
Specializing to the case s = 2 and expanding the integrand, we obtain exactly (2.57) that produces the right answer for the one-point function of the stress-energy tensor.
In (2.88), a Lorentz invariant UV regulator for the sum-integral is implicit (see around (2.33)). More explicitly, assuming that Σ ren (ω n , p) is analytic and does not contain additional poles in ω n , we can follow the derivation of (2.77) to obtain the following contour integral expression at general even spin s, where γ is the same contour as in (2.77) and ϵ p = p 2 + ∆ 2 . Pushing the contour to infinity (we have implicitly suppressed the regulator ϵ in (2.77) and the divergences in ϵ are cancelled by counter-terms), we obtain (2.90) The above expression is manifestly finite but still challenging to evaluate numerically due to the complicated expression for Σ(ω n , p) as in (2.84). 15

Gross-Neveu Model and Variations
We now apply the methods developed in Section 2 to study thermal observables in d = 3 fermionic CFTs with vector-like large N limits. We will focus on the Gross-Neveu (GN) model as well as its closely related variations.

Review of the Gross-Neveu Model and Large N Expansion
We start by reviewing the celebrated Gross-Neveu (GN) model [53] (see also [18]), keeping spacetime dimension d general for the moment. This model describes N interacting Dirac fermions ψ i governed by the following action which has an obvious U (N ) global symmetry that rotates the Dirac fermions, as well as a Z 2 parity symmetry that acts as, 16 which forbids the U (N ) invariant Dirac mass term mψ i ψ i . For d > 2, the quartic interaction is non-renormalizable and this theory needs a proper UV completion. This is provided for 2 < d < 4 [54] by the Gross-Neveu-Yukawa (GNY) model which has in addition a real pseudoscalar φ (which is Z 2 parity odd) and the following action, symmetric and parity invariant critical point is commonly referred to as the GN or the GNY CFT. It describes the second-order order-disorder phase transition with order parameter 15 One may hope to derive a simpler expression for the one-point function coefficient b 1 s as we have done for the s = 2 case in (2.64) (see also (2.40)) directly in terms of the leading self-energy Π β (ω n , p) for the σ field whose explicit form is given in (2.37). 16 This comes from a time-reversal symmetry in the Lorentzian signature after Wick rotation. φ ∼ψψ and characterized by the spontaneous Z 2 parity symmetry breaking and dynamical mass generation for the fermions (see [18] for a more extensive review).
In d = 3, which will be the focus here, the GN ( 1)) ⋊ Z C 2 symmetry (also known as the Nambu-Jona-Lasino model) which we will come to in Section 3.5. They describe various quantum phase transitions for spinless and spinful fermions on lattices including graphene [62][63][64][65]. There has been recent progress in determining the flat space CFT data in these fermionic CFTs using the bootstrap method (see [66] and references therein).
The large N solution of the GN model at criticality can be deduced in a similar way as for the critical vector model. We introduce an auxiliary scalar field ϕ, analogous to the σ field for the O(N ) vector model and implement a Hubbard-Stratonovich (HS)transformation We will set g 0 = N g t 0 to obtain the proper large N limit with g t 0 held fixed. Integrating out the fermions, we arrive at the following effective Lagrangian for ϕ, As explained in [54], the GN and the GNY models coincide in the scaling region where the HS field ϕ coincides with the Yukawa pseudoscalar φ up to a normalization factor. 17 A small subtlety is that for N = 1 2 , the CFT is properly defined only in the GNY description as the four-fermion coupling vanishes for a single Majorana fermion. Interestingly, this model has emergent N = 1 supersymmetry at the fixed point [55] and is also known as the N = 1 super-Ising CFT [56].
In the large N limit we can again argue that the path integral over ϕ is dominated by its saddle-point, which satisfies the so-called gap equation where c d ≡ tr 1 counts the components of the Dirac spinor in d spacetime dimensions (i.e. c 3 = 2 in the case of interest). This equation again contains divergences, that we cancel by introducing a renormalized coupling constant g t , where K d is scheme-independent and given by (2.11) and in the last equality above we have neglected further corrections that are scheme-dependent (and suppressed near the critical point). To reach criticality, we fine-tune g t such that ϕ = 0 so that the fermions become massless. One can check that at this point the HS field ϕ is also massless. Indeed, the inverse propagator for ϕ is, 18 in the dimensional regularization scheme, The second term on the RHS above vanishes as p → 0 and we have G −1 ϕ (0) = −g t . Therefore ϕ is massive unless g t = 0.

Free Energy at the Subleading Order and Chemical Potential Dependence
The thermal free energy density or equivalently the one-point function of the stress-energy tensor in the GN CFT has the same structure as a series in 1 N as in (2.21) for the critical O(N ) model. In this section we will compute the leading contribution F GN,0 ∝ N and the first subleading correction F GN,−1 ∝ N 0 in this expansion. The computation of both contributions is similar to the critical O(N ) model. We will comment on the differences and the subtleties associated with these fermionic CFTs along the way.
We start with the leading contribution F GN,0 (β). We assume that at finite temperature, the dominant configuration for the HS field ϕ is homogeneous and thus the free energy density is where the saddle point value of ϕ will be fixed shortly. There are two spin-structures for the fermions that are compatible with the O (2N ) global symmetry, corresponding to either periodic or anti-periodic boundary conditions along S 1 β . Correspondingly the allowed frequencies in the sum of (3.10) are, periodic : ω + n = 2πn β , antiperiodic : ω − n = 2π(n + 1/2) β , (3.11) and the latter is the standard thermal boundary condition for fermions. In the periodic case, the free energy density (3.10) as a function of ϕ can be read off from (2.28) by comparing (3.10) and (2.22), One finds from (2.31) that the saddle point and the final free energy density are Note that F P GN,0 is positive, naively suggesting that the specific heat of such model is negative, but we are considering the periodic boundary condition, which does not define a conventional thermal ensemble.
In the anti-periodic case, following a similar derivation to (2.28), the free energy density before plugging in the saddle point for ϕ reads (see also [10,11] (3.14) Formally, F AP GN,0 (ϕ) has three critical points, The last two solutions ϕ = ±ϕ * are equivalent by the Z 2 parity symmetry and have lower free energy density than the first solution. Furthermore, by considering subsequent 1 N corrections, we can check that ϕ = 0 is incompatible with conformal symmetry. Therefore to describe the large N thermal CFT, we pick ϕ = ϕ * as our saddle-point solution. Then the leading contribution to the free energy density is where the Clausen function is related to usual polylogarithms by We now determine the subleading 1 N corrections to the free energy density in the GN CFT, by taking into account the fluctuations of ϕ around its saddle point. This is related to the computation of the propagator for ϕ in such a non-trivial background. In this case, as shown in Appendix C, the ϕ propagator takes the following form, where P = (Ω, p) is a short-hand for the three-momentum with P 2 = Ω 2 + p 2 , and ϕ + ≡ √ σ * , ϕ − ≡ ϕ * are the saddle points. For fermions obeying the periodic boundary condition Π + β (Ω, p) = Π β (Ω, p) as is given in (2.36) (see also (2.37)) and for the anti-periodic case Π − β (Ω, p) is given in (C.6). The subsequent 1 N correction to the free energy density of the GN CFT with the periodic boundary condition is related to the O(N ) answer as follows, where we have already computed the second term in the last line in (2.46) and the first term after proper renormalization follows from (2.28), We thus arrive at the final answer of the subleading free energy density for the GN CFT in the periodic spin structure, For the anti-periodic spin structure, an analogous computation yields the following expression for the self-energy factor for ϕ (for derivation see Appendix C.2), , (3.22) and the contribution to the free energy density is At first glance, for ϕ 2 * < 0 the integrand becomes imaginary at some range of energymomentum, and therefore seem unphysical. However, we should remember that the contour of integration for the field ϕ needs to be deformed (to be the steepest descent contour) in order to implement a proper saddle-point approximation. To take that into account, we just need to take the absolute value of the argument of logarithm (that corresponds to the slight tilt of the contour of integration), we then arrive at where we have used from (2.28), 25) and also (2.29). Note that it is very important that we have chosen the right solution among and introduce an external constant electromagnetic potential A 0 = µi β with µ ∈ [0, 2π). It amounts to shifting frequencies of the fermions to βω n ≡ βω − n + µ. Note that in particular µ = π corresponds to the periodic spin structure for the fermions. The leading free energy density as a function of µ and the one-point function of ϕ field (denoted byφ below) is Then one can compute the propagator of ϕ field on such a thermal background with chemical potential (see (C.7) for the explicit expression) and then determine the subleading correction to the free energy density F GN,−1 (β, µ) as a function of chemical potential µ (see Figure 6).

One-point Functions of Higher-Spin Currents
As explained before, the free energy density of the CFT determines the thermal one-point function of the stress-energy tensor. The stress-energy tensor belongs to a tower of spinning operators known as the higher-spin currents, which in the GN model are given by the following, with s ≥ 1, where e µ is a unit vector along the S 1 β and∆ s is the scaling dimension of the (weakly broken) higher-spin current,∆ s = s + 1 +γ s , (3.30) with the anomalous dimensionγ s that is suppressed by 1 N and can be found in [52,[68][69][70]. The ± subscript in (3.29) refers to the choice of spin structure along the S 1 β . To compute the coefficientsb ±,s in the leading large N limit, we introduce an auxiliary null three-vector ξ as in Section 2.6 (see around (2.69)). Following the same procedure there with the same type of Feynman diagrams (the propagating bosons are substituted by the fermions), we find that where for ϕ + = √ σ * we sum over βω + n = 2πn, and for ϕ − = ϕ * we sum over βω − n = π(2n + 1) with n ∈ Z. The UV divergences in the above expressions are regularized as for the critical O(N ) model. For the anti-periodic boundary condition, we consider the following generating function similar to (2.72), Expanding the RHS of the above equation we find thatb 0 −,s = 0 for s odd as expected and for s even we have,b We can also derive an alternative formula forb 0 −,s starting from (3.31). Following the procedure around (2.77), we arrive at the following integral expression (setting β = 1 below), for where ϵ p ≡ p 2 + ϕ 2 * . The above integral reproduces the sum expression (3.34) thanks to the following integral identity, The subleading 1 N correctionsb 1 ±,s to the one-point functions of these higher-spin currents in the fermionic CFT can be in principle computed in a similar way as describe near the end of Section 2.6 for the O(N ) scalar model.

The Chiral Ising Gross-Neveu-Yukawa Model
As it was explained in the previous section, another interesting generalization of the usual GN or GNY model is the chiral Ising GN or GNY model with a different symmetry, which has many applications in condensed matter systems including the semimetal-insulator transition in graphene [62][63][64]. The chiral Ising GN model is defined by splitting the N Dirac fermions into two groups as ψ i = (χ L a , χ R a ) where a = 1, . . . , N 2 with the following action that has a different quartic interaction compared to (3.1), 22 Consequently the global symmetry of the model is reduced where the so-called "chiral" Z chiral 2 exchanges χ L a and χ R a . In d = 3, this symmetry is enhanced to O(N ) 2 ⋊ Z chiral 2 as mentioned previously.
For 2 < d < 4, the UV completion of the model (3.37) is provided by the chiral Ising GNY model, with one additional real pseudoscalar φ and the following action, symmetry at the semimetal-insulator transition in graphene for which φ is the order parameter [62][63][64] (see also [66]). In the scaling region, the chiral Ising GN and GNY models are expected to coincide [54].
Implementing the HS transformation, the action (3.37) can be recast into the following where the HS field ϕ is again identified up to an normalization with the pseudoscalar φ in the chiral Ising GNY model (3.38). We then integrate out the fermions and arrive at the 22 Again there is little difference to the description of the model when N is odd for spacetime dimension d = 3, in which case we simply write the action in terms of the Majorana fermions.
following effective Lagrangian for the field ϕ F cGN (ϕ) = −N tr log (i/ ∂ + ϕ) + tr log (i/ ∂ − ϕ) where g 0 = g t 0 N as before with g t 0 held fixed in the large N limit. Note the two terms in (3.40) where ϕ comes with opposite signs. For a general configuration of ϕ, these two terms are not related to each other and that's why the chiral Ising GNY (or GN) model is different from the usual GNY model. Nonetheless, as we explain below, it is easy to show that the leading and the first subleading corrections to the free energy density coincide. 23 This is because for constant ϕ, the first two terms in (3.40) are identical and the saddle point of ϕ obeys the same gap equation as for the GN model, and by fine-tuning the coupling g t 0 we can bring system to criticality as before. The free energy density of the chiral Ising GN model has the same large N expansion as before, where F cGN,α scales as N 1+α . The above discussion indicates that the leading N contributions agree F ± cGN,0 = F ± GN,0 . Furthermore, the first subleading corrections also coincide F ± cGN,−1 = F ± GN,−1 because the large N propagator of the field ϕ is an even function of its constant saddle point value (see (3.18)). However it is easy to see that already at the second order, there is a difference between the two versions of the GN model in the free energy density, namely F ± cGN,−2 ̸ = F ± GN,−2 , which arises from the different contributions of the vacuum diagrams in Figure 7 at non-zero temperature. The same conclusion holds for the thermal one-point function of general higher-spin currents defined in (3.28).

The Nambu-Jona-Lasinio (NJL) Model
Another important generalization of the GN model reviewed in Section 3.1 is the Nambu-Jona-Lasinio (NJL) model. It is commonly defined in terms of four-component Dirac fermions 23 In fact the same analysis below easily extends to the free energy on S 1 β × S 2 R and the same conclusion that the free energies for the GN model and its chiral version only differ at the 1 N order continue to hold. It would be interesting to study this difference in free energies in relation to the difference in the operator spectrum as was analyzed in [66]. Figure 7: The first two diagrams in the 1 N expansion whose contributions to the free energy density differ between the GNY and the chiral Ising GNY models (due to the choice of signs at the cubic vertices in the latter). The solid lines denote the propagators for the fermions ψ i and the dotted lines for the pseudoscalar φ.
Ψ i with i = 1, . . . , N f (which naturally arise in the d = 4 model) with the following action, where γ 5 is the usual chirality matrix for d = 4 Dirac spinors. 24 The global symmetry is the obvious flavor symmetry that rotates the Dirac fermions Ψ i , the extra U (1) "chiral" symmetry acts by Ψ i → e iαγ 5 Ψ i , and Z C 2 is the charge conjugation symmetry. For d = 3, the global symmetry is enhanced to As is the case for the other fermionic models with four-fermion interactions, the model (3.43) is non-renormalizable for d > 2 and its UV completion for 2 < d < 4 is provided by the Nambu-Jona-Lasinio-Yukawa (NJLY) model which contains one complex scalar field φ = φ 1 + iφ 2 and has the following action [54], The scalar field transforms under the chiral U (1) symmetry as φ → e iα φ. In the scaling region, the NJL and the NJLY models are expected to coincide [54], where the theory undergoes a second-order phase transition with order parameter φ and spontaneous U (1) chiral symmetry breaking.
The NJL (NJLY) model is also known as the chiral XY GN (GNY) model in the condensed matter literature, where it describes various quantum phase transitions with complex order parameters (see for example [72][73][74][75][76]). In particular, it governs the superconducting phase transition in graphene at N f = 2 [72,73] and the semimetal-VBS (valence bond solid) transition in graphene and graphene-like material for various values of N f [74] as well as the critical surface states of certain topological insulators with emergent supersymmetry at [73,77,78]. 25 To solve the NJL model in the large N = 2N f expansion, we proceed by introducing a complex HS field ϕ = ϕ 1 + iϕ 2 , which is identified (up to an overall normalization) with the complex scalar φ in the NJLY model near the critical point. We then integrate out the Dirac fermions and obtain the following effective Lagrangian, with where the rescaled coupling g t 0 is related to the coupling constant in (3.43) by g 0 = N g t 0 and kept fixed in the large N limit. Again by fine-tuning g t 0 we can bring the model to its critical point which is described by the chiral XY GN (or GNY) CFT. We exploit the semi-classical approximation in the large N limit as before and derive the gap equation that governs the dominant saddle which we assume to be described by a homogeneous configuration of the complex HS scalar ϕ. Taking advantage of the U (1) symmetry, we can always set ϕ 2 = 0. At zero temperature, the gap equation is then equivalent to the one for the ordinary GN model Nonetheless, as we show explicitly below, they start to differ at the first subleading order in 1 N . To obtain the 1 N correction F P,AP NJL,−1 to the free energy density, we compute the induced propagators for the HS fields ϕ 1 , ϕ 2 , which take the following form, with ϕ ± and Π ± β (Ω, p) the same as around (3.18) where the ± label the periodic and anti-periodic spin structures respectively. Note that the additional zero of G ± ϕ −1 22 at P 2 = 0 corresponds to the Goldstone boson for the chiral U (1) symmetry, which contributes to the 1 N corrections to the free energy density of the critical NJL model. Comparing with the analysis in Section 3.2 for the GN model, we find that the full 1 N correction in the NJL model for the free energy density is given by, for the periodic spin structure, and analogously for the anti-periodic (thermal) spin structure, the result is As promised, the above clearly differ from the results in the GN model (see (3.21) and (3.24)).
Let us know consider more general thermal observables given by one-point function of the higher-spin currents in the critical NJL model. They are defined in a way analogous to (3.28) but in terms of the four-component Dirac fermions, The corresponding one-point function coefficientsb ±,s (defined in a similar way as in (3.29)) can be computed from Feynman diagrams as in Section 3.3. As one may expect from the discussion around (3.46) that concerns the free energy density (equivalently the stress-energy tensor one-point function), these more general one-point functions in the critical NJL model also coincide with those in the ordinary GN model (see (3.32) and (3.34)) to the leading order in the 1 N expansion, namelyb 0 ±,s =b 0 ±,s .
We expect this agreement to fail at the subleading order which we have demonstrated explicitly for s = 2 (see around (3.48) and (3.49)) and for general s by an analysis parallel to that in Section 2.6.

Chern-Simons Quantum Electrodynamics
We now study finite temperature observables in d = 3 vector-like large N CFTs whose definition involves dynamical gauge fields. For concreteness, we focus on perhaps the simplest class of nontrivial gauge theories, namely the quantum electrodynamics (QED) with N massless Dirac fermions and Chern-Simons coupling k (which we refer to as CS QED), and present explicit results for its free energy density at the subleading order in the 1 N expansion as a function of the 't Hooft coupling λ = 4πN k .

Review of the CS QED and Large N Thermal Observables
The CS QED model in d = 3 describes the interaction of N two-component Dirac fermions with the electromagnetic field A µ , with Maxwell coupling e 0 , fermion bare mass m 0 and Chern-Simons coupling k that satisfies the quantization condition k + N 2 ∈ Z. 26 The theory has a U (N )/Γ ⋊ Z C 2 global symmetry where the U (1) subgroup of U (N ) comes from the topological symmetry (under which the monopole operators are charged), the SU (N ) subgroup is the flavor symmetry of the Dirac fermions, Z C 2 is the charge conjugation symmetry and the discrete quotient Γ ⊂ U (1) ensures that the symmetry acts faithfully on the local operators [81,82]. 27 There is substantial evidence that at sufficiently large N , upon tuning the bare mass m 0 , the 3d CS QED flows to a nontrivial U (N ) invariant CFT at long distance [83][84][85][86][87][88], which we will refer to as the conformal CS QED. 28 This is particularly well-established by considering the limit of N → ∞ and studying the 1 N expansion (see for example [12,[95][96][97][98][99]). 29 See 26 Note that the mass term in (4.1) is parity-odd and SU (N ) invariant. The CS level k here (also known as the effective CS level) includes the formal N 2 shift from the gauge invariant regularization of the fermion determinant using the eta invariant exp N 2 πiη(A) [80]. The theory is parity invariant only when N is even and the effective CS level vanishes k = 0. 27 For k < N 2 , Γ = Z N/2−k and otherwise Γ = Z N/2+k . See [81] for details. 28 These abelian gauge theories are also related by conjectured IR dualities to certain nonabelian Chern-Simons-matter theories [89][90][91][92][93][94]. Consequently our results for the abelian theories in the large N limit give rise to predictions for the corresponding thermal observables in these nonabelian theories, which are strongly coupled, due to the very nature of these level/rank type dualities. 29 At small enough even N ≤ N * for even critical flavor number N * , the QED model with vanishing Chern-Simons level k = 0 is believed to exhibit spontaneously "chiral" symmetry breaking of SU (N ) → SU (N/2) × SU (N/2) × U (1) due to certain parity-even SU (N ) invariant scalar operator (quartic in the fermions) becoming relevant and inducing a nonzero expectation value for the order parameter (quadratic in the fermions) [83,84,88,100]. also [101][102][103][104] for recent works on bootstrapping the corresponding CFTs directly at finite N . Such CFTs are known to describe exotic phases of quantum matter in two spatial dimensions including the Dirac spin liquid [105][106][107][108][109][110] for vanishing Chern-Simons level k = 0 and various values of N depending on the underlying two-dimensional lattice. For nonzero Chern-Simons level k ̸ = 0, they also emerge near the quantum phase transitions between different quantum hall states [111]. Below we explain how to extract finite temperature observables for these quantum critical points in the large N limit, taking into account the subleading effects.
By integrating out the fermions, we arrive at the following effective Lagrangian for the gauge field, where we have rescaled the couplings and we take the large N limit with g 0 and λ kept fixed. Note that the gauge field A µ here plays the role of the ϕ field in the case of the GN model (similarly the σ field for the scalar It is easy to see in the large N limit that the theory has a nontrivial fixed point in the IR for any initial coupling g 2 0 . The quantum correction to the fermion mass in the IR can be cancelled by properly fine-tuning the bare mass parameter m 0 (and in the dimensional regularization scheme this corresponds to setting m 0 = 0). 30 The propagator for the gauge field A µ can be computed from standard Feynman diagrams (see Figure 8) and in the Lorentz gauge is given by (see also [12]), and at large distance, which indeed takes the right form that is required for conformal symmetry. In particular, 30 No fine-tuning is required for k = 0 and N even due to the extra parity symmetry.
it implies that the gauge invariant operator F µν is a primary operator of scaling dimension ∆ = 2 in the large N limit. This also means that in the IR we can neglect the Maxwell action since it is irrelevant. We proceed as before to study this CFT on S 1 β × R 2 to determine its thermal observables. The free energy density of the conformal CS QED has the same type of 1 N expansion as in (2.21) for the critical O(N ) model, and we will be interested in the first two terms F CSQED,0 (β) at order O(N ) and F CSQED,−1 (β) at order O(N 0 ).
We can again argue that in the large N limit we can pursue a semi-classical approximation where A µ is homogeneous. Using gauge transformations we can set the spatial components along R 2 to A 1 = A 2 = 0, but we can not set A 0 = 0 because of the non-trivial topology of the thermal background. Nonetheless, we can fix and compute the effective Lagrangian as a function of the holonomy u ∈ [0, 2π/β), where we have imposed anti-periodic boundary conditions 31 for the fermions along S 1 β . Implementing the regularization for the sum-integral as before (see in particular (3.26)), we arrive at, (4.8) The homogeneous saddle point for the holonomy u relevant for the thermal CFT and the 31 This is without loss of generality since u changes the effective boundary condition for the fermions.
corresponding leading free energy density then follow, which coincides with the free energy density of N free Dirac fermions. We emphasize that F CSQED (u) has another saddle point at u = π which is a local maximum and thus unstable.
We see that effective anti-periodic boundary conditions on the fermions are energetically favored.
It is easy to repeat the analysis in the previous section to compute the one-point function of the higher-spin currents defined in (3.28) in the large N limit of the conformal CS QED.
The one-point function coefficientb 0 s for the spin s current follows from (3.34) with ϕ * = 0, which again coincides with the free fermions as expected in the leading large N limit of the CS QED.

Subleading corrections to the Free Energy
The large N CS QED is of course very different from free fermions. Indeed, the distinctions already manifest at the first subleading order in the 1 N expansion of thermal observables in the CFT. Below we demonstrate this explicitly for the free energy density (equivalently the one-point function of the stress-energy tensor), which depends nontrivially on the 't Hooft coupling λ.
To obtain the 1 N correction to the free energy density, we compute the quantum propagator for the gauge field A µ . We first work with vanishing Chern-Simons level k = 0. The photon polarization tensor then reads, , Ω n ≡ 2πn β , (4.11) from fermions running in the loop. One can check that it has the following structure as required by the residual symmetry at finite temperature and the Ward identity, Π 00 (Ω n , q) = Π E (Ω n , q) , Π 0i (Ω n , q) = − Ω n q i q 2 Π E (Ω n , q) , where Π E (Ω n , q) and Π M (Ω n , q) are scalar functions (see Appendix D for explicit expressions). After taking into account the Faddeev-Popov ghosts, the 1 N correction to the free energy density of the conformal QED is given by, (4.13) Using the results of the Appendix D and implementing the numerical evaluation of the sum-integral above as explained in Section 2.3, we find that, which agrees with the results in [12]. We now re-introduce the Chern-Simons coupling k and study the dependence of the free energy density on the 't-Hooft coupling constant λ = 4πN k . Recall that the leading contribution to the free energy density F CSQED,0 (β) is λ independent. In this case, the photon polarization operator is modified in a simple way to, Consequently, the free energy density at the subleading order reads, which now depends nontrivially on λ and can be computed numerically from the method outlined in Section 2.3 (see Figure 9). In particular, it interpolates between the free fermion answer at λ = 0 where there is no subleading correction to the free energy density (4.9) and the pure QED answer (4.14) at λ = ∞.

A Lattice Regularization
Here we discuss an alternative approach to compute the thermal free energy of the critical O(N ) model in d = 3 using a finite lattice. First, we consider the system on an infinite lattice Z 3 with lattice spacing a [112]. It amounts to replacing the scalar propagators in the following way with normalized periodic momenta p x,y,z ∈ [−π, π]. This automatically regulates the UV divergences. The gap equation (2.6) becomes, expanding the left-hand side of this equation for small σ we obtain, where λ t * is determined from numerically fitting the integral. To approach the critical point we fine-tune r t 0 to r t * . In this case σ = 0 is the solution and the system is gapless in the IR. The other coupling λ t 0 could be left arbitrary but it is better for the numerics to also tune it to λ t 0 = λ t * to cancel the leading finite-size corrections from the least irrelevant operator. Indeed, the propagator for the HS field σ in this limit becomes, and if we set λ t 0 = λ t * the leading correction to the conformal propagator is cancelled. Now we study the same model at criticality on a finite lattice Z nt × Z 2 n . For general n, n t ≫ 1 this correspond to the lattice approximation for studying the CFT on the torus If we further take n ≫ n t , the lattice setup approximates the thermal background R 2 × S 1 β with inverse temperature β = n t a. The scalar propagator on this lattice is G −1 a (p, σ) = 6 − 2 cos p x − 2 cos p y − 2 cos p t + σ , p x,y = 2πi x,y n , p t = 2πi t n t .
(A. 6) where i x,y,t ∈ Z. The gap equation is and we will attempt to set λ t 0 ∼ λ t * to reduce the leading finite-size corrections as discussed above. We have checked numerically that the solution of this gap equation is σ * = ∆ 2 n 2 t which produces the leading free energy density (see (2.32)) with good precision for a lattice of size n = 100 and n t = 40 To obtain the subleading correction to the free energy density, we should study the propagator of the HS field σ, equivalently the self-energy (see also (2.36)), (A.9) In the limit of the interest n ≫ n t ≫ 1 and small q, we expect that Π β (q) ∼ n tΠσ (qn t ).
Then the subleading correction to the free energy is Limited by the precision for the value of λ t * in (A.4) and the size of the lattice we can consider with the available computing power, we wouldn't get the desired behavior for the free energy density at the critical point (where f −1 is the desired CFT observable and the constant needs to be subtracted by numerical fitting).
To see this more explicitly, we estimate that the self-energy Π β (q) behaves for n t ≫ qn t ≫ 1 as (A.12) To achieve a good approximation to the critical behavior, we need to consider a sizable lattice such as n ∼ 500 , n t ∼ 25. For instance, in this parameter regime, the value of the σ field from solving the gap equation is σ * n 2 t = 0.9268 which is close to the critical value (2.31). But then to evaluate the sum in (A.10) accurately requires large processing powers that we don't have at the moment. For a smaller lattice such as n ∼ 70 and n t ∼ 10 however, we obtain instead σ * n 2 t = 0.9409, which reflects a relative error of ϵ σ * ∼ 1%. Since the free energy density involves a sum over a logarithm of the self-energy which depends on σ, such a relative error is quite important.
For smaller lattices, the subleading corrections in (A.12) also need to be taken into account, which requires setting λ t 0 = λ t * with high precision. However the fact that we don't know λ t * (see (A.4)) at high precision introduces a numerical uncertainty to the computation of the Π β which again propagates to the numerical evaluation of the subleading piece F −1 (n t ) of the critical free energy density. We expect the total error in (A.11) to be ∆ F −1 ∼ 0.01 for n t ∼ 10, which makes it hard to extract the small critical amplitude f −1 (see (2.46)). To improve the situation, we must know λ t * at least within a relative error of ϵ λ t * ∼ n −3 t .

B Self-energy of σ in the O(N ) Vector Model
In this section we analyze the two-point function of the HS field σ at finite temperature for the critical O(N ) model at large N and provide useful analytic expressions that will help the numerical evaluation of subleading effects for thermal observables in the CFT.
From diagrammatic analysis (equivalently by saddle-point approximation of the path integral), one finds that the self-energy for the σ field (i.e. the inverse propagator) is, where ω n ≡ 2πn β denotes the Matsubara frequency and σ * is the solution to the gap equation given in (2.31). Using Feynman parametrization, we take the integral over spatial momentum and obtain the following, where we have introduced M 2 ≡ σ * + (Ω 2 + q 2 ) x(1 − x) > 0. Evaluating the sum in (B.2) then gives, x)) .

(B.3)
If we introduce a chemical potential µ for the U (1) ⊂ O(N ) subgroup with even N as described in Section 2.4, the self-energy becomes, up to regularization and renormalization as described in Section 2.2.

B.1 Large Momentum Expansion
To analyze the UV structure of the self-energy and the sum-integral in (B.5), it is useful to have the large momentum expansion of (B.1) which we provide here, generalizing the results in [9] to further subleading orders.

C Self-energy of ϕ in the Gross-Neveu Model
In this section we compute the two-point function for the HS field ϕ ∼ψ i ψ i in the critical Gross-Neveu model (equivalently the self-energy of ϕ). We consider both periodic and antiperiodic boundary conditions for the fermions.