The unitary Fermi gas at large charge and large N

We study the unitary Fermi gas in a harmonic trapping potential starting from a microscopic theory in the limit of large charge and large number of fermion flavors N. In this regime, we present an algorithmic procedure for extracting data from perturbation theory, order-by-order, without the need for other assumptions. We perform a gradient expansion in the interior of the particle cloud, sufficiently far from the cloud edge where the particle density drops rapidly to zero. In this latter region we present the first microscopic computation characterizing the contribution of the edge terms. The microscopic theory reproduces the predictions of the superfluid EFT, including the action, the form of the gap equation, and the energy of the system in a harmonic trap (which maps, via the non-relativistic state-operator correspondence, to the scaling dimension of the lowest operator of charge Q). We additionally give the Wilsonian coefficients at leading order in N up to NNLO in the large-charge expansion.


Introduction
The large-charge expansion has been shown to lead to important simplifications in strongly-coupled relativistic conformal field theories (cfts) [1][2][3][4].In cases with a unique vacuum at large charge, the physics is that of a superfluid, and an effective field theory (eft) can be formulated as an expansion in inverse powers of the charge.This eft can be used for example to compute the energy of the ground state on the cylinder R × S d in a sector of fixed charge Q, which gives, via the state-operator correspondence, the scaling dimension of the lowest operator of charge Q.
The large-charge expansion also renders useful simplifications for nonrelativistic conformal field theory (nrcft) [5,6].Here, the underlying symmetry is not the conformal algebra but the Schrödinger algebra (summarized, e.g., in [5][6][7]), which naturally contains a conserved U(1) charge (namely, particle number).As in the relativistic case, the symmetry completely determines the form of the two-point function up to the scaling dimension.Moreover, the notion of a state-operator correspondence persists for nrcft: The scaling dimension of a given operator is the same as the energy of the corresponding state for the system in a harmonic trapping potential [7].
The unitary Fermi gas is an interesting and experimentally accessible example of a system with Schrödinger symmetry.It can be realized in the laboratory via an ultracold Fermi gas in an optical trap [8,9].By tuning an external magnetic field to a Feshbach resonance, the scattering length becomes infinite at the so-called unitary point, where the system acquires Schrödinger symmetry.Ultra-cold Fermi gases in the presence of a broad Feshbach resonance can be described by a contact interaction.One requires that the gas is sufficiently dilute, which is the typical case in experiments [10].This amounts to the inter-particle distance being much larger than the range of the inter-atomic potential, so that the details of the interactions are not important, and it is safe to focus on s-wave scattering. 1 Crucially, the optical trap can be chosen to provide the required harmonic trapping potential for the realization of the state-operator correspondence.Thus, we have an accessible laboratory system that, insofar as finite-charge effects are small, matches our theoretical set-up and invites direct comparison of our predictions to experimental data.One example of a measurable observable is the doubly-integrated axial density, which was determined in [11].
More in detail, the unitary Fermi gas can be described by a non-relativistic superfluid.The structure of the effective action was first worked out by Son and Wingate [12], up to a set of undetermined Wilsonian coefficients.The same eft was later realized through the technology of the large-charge expansion [5,6,13,14].To next-to-leading order (nlo) 2 , 1.Some nuclear systems feature a BCS-BEC crossover similar to that of cold Fermi gases, which can be described by the same model, though the interaction must have a finite range, and the crossover is realized by varying the density, as opposed to the scattering length [8]. 2. Sub-leading corrections to the leading-order effective Lagrangian can come from two channels: Derivative in 3 + 1 dimensions, the effective Lagrangian with external potential V(r), takes the form [12] where θ is the Goldstone of the U(1) superfluid, and When evaluated on the classical ground state θ = ℏµt, we obtain As an aside, this form of the effective action can be obtained purely from locality and dimensional arguments.The charge density on the ground state is given by . (1.4) In the absence of a potential, only the first term in ℒ GS survives, and the behavior of the Fermi gas is captured by the single Wilsonian coefficient c 0 .Traditionally, this coefficient is expressed via the so-called Bertsch parameter ξ [12,15], defined as the ratio of the energy density of the unitary Fermi gas to the energy density of a free Fermi gas at the same density, so that Once a potential is included, new scales appear and the physics is no longer captured by a single Wilsonian coefficient (equivalently, the Bertsch parameter), as defined above.
For the harmonic trap, the bulk contribution to the ground-state energy, respectively the scaling dimension ∆ of the lowest operator with charge Q, is given by where Q is the fixed particle number and N is the number of fermion flavors.When considering the set-up of a cloud of particles in a harmonic trapping potential, corrections and loop corrections.The former contribute to what are labeled as NLO terms in (1.1), while the latter are further suppressed.Thus, to the order of interest, tree-level analysis will suffice.
it is obvious that the above eft is only valid in the bulk, sufficiently far away from the cloud edge located at a distance R cl from the center, where the particle density rapidly drops to zero.For this reason, the bulk eft has to be complemented with terms from the edge [16,17].The most general form of terms located at the cloud edge for a system in d + 1 dimensions is where p is an integer, κ p is a Wilsonian coefficient, and δ(U) is an operator-valued delta function that localizes on the cloud edge.Generally speaking, the contributions of the bulk operators to ∆ can have divergences stemming from the cloud edge if d is even and the operator has positive Q-scaling.These divergences can always be regulated, however, by an edge counterterm of the same µ-scaling, leading to log(Q) terms in ∆ [14,16,17].Additionally, there is a universal log(Q) term in odd d descending from the Casimir energy of the Goldstone θ.
The eft result for ∆, including contributions from the bulk, the edge, and the Casimir energy, specialized to 3 + 1 dimensions, is given by [14,17] Here, the terms starting with Q 12/9 come from bulk operators, those starting with Q 5/9 come from edge operators, and log Q is a universal contribution from the Casimir energy.Further terms are suppressed by negative powers of Q.The coefficients a i , b i are related to the Wilsonian coefficients and cannot be computed within the framework of the eft alone.
In this work, we follow a different approach.We start from a microscopic action, including the trapping potential, and find the form of the effective action in the limit of large N fermion flavors, in a sector of large charge and zero temperature.As required by locality and the symmetries of the system, we indeed reproduce the form of the ground-state effective action in Eq. (1.3).Importantly, however, the specific microscopic action allows for the direct computation of the Wilsonian coefficients c 0 , c 1 and c 2 of the bulk eft.
In principle, can we follow the same strategy as in the relativistic case of the O(N) model at large charge and large N [18], where standard large-N path integral techniques such as the Stratonovich transform can be employed in a sector of fixed charge.Due to the inhomogeneity introduced by the potential, however, we require more sophisticated techniques.In fact, evaluating the functional determinant from the Gaussian integral of the original fermionic degrees of freedom (dof) becomes a nontrivial problem in the presence of the potential.Only when working at large charge is it possible to perform a gradient expansion in the bulk, as the dimensionless parameter ε ∼ 1/µR 2 cl in front of the gradient plays the role of ℏ in the quantum mechanical problem of the heat kernel.At large charge, there is a spontaneous symmetry breaking (ssb), which, at large N, is due to an explicit realization of the Cooper mechanism.This is because the collective field σ(x) (which acquires a non-trivial vacuum expectation value (vev) and plays the role of a gap for the fermions) is a bilinear of the fundamental fermions, (1.10) Below, we compute the form of the gap to second order in 1/ε, finding with y 0 ≈ 1.1622 . . ., y 1 ≈ −0.00434691 . . ., and y 2 ≈ −0.160794 . . . .We also compute the energy in the harmonic trap to second order in 1/ε, (again, corresponding to the scaling dimension of the lowest operator of charge Q), From there, we can extract the Wilsonian coefficients in the bulk, c 0 , c 1 , c 2 , c 0 = 0.0841 . . ., c 1 = 0.006577 . . ., c 2 ≈ 0.004872 . . ., (1.13) and the value of the Bertsch parameter, ξ ≈ 0.5906 . . . .(1.14)This bulk expansion, however, is no longer valid once we approach the cloud edge.Unfortunately, at the edge, we have no small parameter in which to expand the functional determinant (see Figure 1).We are, however, able to estimate the scaling of the edge contributions, which indeed matches the eft prediction [16,17].
In the previous literature, a number of calculations have been made using various different approximations and approaches.In [19], a large-N computation was made without an external potential, including subleading corrections in 1/N.Our leadingorder result (the Bertsch parameter) agrees with them.A similar computation has been done in [20], which also addresses the case of finite temperature.A mean-field study was performed in [21].A direct comparison of our results can be made with the work of [22], which gives the corrections to gap equation and the charge density in 1/ε 2 .We are in complete agreement for the shape of the gap.There appears to be a discrepancy at next-to-next-to-leading order (nnlo) between our Wilsonian coefficients and those extracted from their expression for the charge density, however this difference is due to a total derivative term.There are also numerical lattice [23] studies, though these deal with the opposite regime of small charge, where the quantized nature of the problem is apparent, and it is not clear how to make a direct comparison between our respective results.
The main difference with these previous works is that we use a well-justified approximation, within a controlled perturbation theory, and whose only constraints amount to standard dimensional and symmetry arguments.This renders a clean, algorithmic way of extracting perturbative data from a very basic set of assumptions.(In fact, the large-charge, large-N regime justifies the commonly used Thomas-Fermi approximation.)Moreover, we are able to extract the parametric form of the edge terms from a microscopic description.To our knowledge, this is the first time the droplet edge has been tackled from first principles.
The plan of this article is as follows.In Section 2, we discuss the microscopic model, which is our starting point at large N, and use standard large-N techniques to compute a functional determinant upon performing the Stratonovich transform and integrating out the original fermionic dof.In Section 3, we introduce Wigner coordinates and the Moyal product, which allows us to compute the large-charge expansion of the heat kernel in the bulk up to nnlo.In Section 4, we discuss the edge expansion.In Section 5, we give conclusions and an outlook.In the appendices we provide an alternative regularization for the grand potential.

The model
Building upon [19,20,24,25], we consider a model of N fermion flavors with two hyperfine ("spin") states ψ σa , where a ∈ {1, . . ., N} and σ ∈ {↑, ↓}, with a unique chemical potential µ, coupled to an external potential V(r).Working in Euclidean time τ at zero temperature, i.e., on R 4 , this system is described by the action where ψσa , ψ σa are Grassmann fields 3 and Summation over repeating flavour and spin indices is assumed throughout.We keep the potential generic, as we are interested in computing the Wilsonian coefficients in the eft. 4  The bare coupling u 0 < 0 describes an attractive 4-fermion contact interaction which, after renormalization, is proportional to the s-wave scattering length a s .The latter is an experimentally tunable parameter, which diverges at the so-called Feshbach resonance.This point in the phase diagram of the Fermi gas is called the unitary limit, because the s-wave scattering cross-section saturates the bound imposed by unitarity of the S-matrix [8].In this limit, the sole physical length scale characterizing the system drops out.There is strong evidence that this emergent scale invariance is preserved at the quantum level, thereby constituting an example of a strongly interacting nonrelativistic cft.
The action has a manifest (U(1) ⊗ SU(2)) N symmetry, where the U(1) copies correspond to particle number conservation of every flavor.However, the complete symmetry group is larger.Indeed, while the kinetic term has an explicit U(2N) symmetry, the interaction 3. ψσa denotes the complex conjugate of ψ σa .Note that we work in units where k B = 1, so at finite temperature T = β −1 , these fields would have antiperiodic boundary conditions along the thermal circle S 1 ℏβ , with Matsubara frequencies ω n = (2n + 1)π/(βℏ) (n ∈ Z) [26].
4. Special forms of V(r) at this stage could lead to accidental symmetries and identifications of operators in the effective action, for instance.
term written in the form where The true symmetry group of the action is thus U(1) ⊗ Sp(2N), with the usual definition [25].Since the Cooper interaction (2.3) appears explicitly in the bare action, we expect ssb of the global U(1) symmetry to take place at fixed charge.This is to be contrasted with the Gross-Neveu model [27], where the breaking is exponentially suppressed (and thus nonperturbative) at large N.
We can introduce a so-called (Hubbard-)Stratonovich field σ(τ, r) in the Cooper channel (see [28] for a discussion of the different channels), namely, resulting in ) where, in the second line, we have used the Nambu representation and Ψa ≡ Ψ † a , along with the inverse fermion propagator The Stratonovich field σ is a bilinear of the fundamental fermions, which is charged under the U(1).Therefore, if it receives a non-zero vev, the symmetry is spontaneously broken.This is the mechanism underlying Bardeen-Cooper-Schrieffer (bcs) mean-field theory.
With the fermions appearing only quadratically in this formulation, we can integrate them out to find the effective action for σ: As usual in models with a vector symmetry, the fact that the quadratic action for σ is proportional to N means that we can separate σ into a saddle-point value σ 0 plus quantum fluctuations suppressed by 1/ with δS[σ] δσ(τ, r) In bcs theory, the saddle σ 0 is usually referred to as the gap (especially in the case without a trap, where it is homogeneous), and the above condition is dubbed the gap equation.In practice, we shall solve it to find the σ 0 -profile only at the very end, therefore treating σ 0 as a generic function.
In Eq. (2.8), the leading-order term in a 1/N expansion, namely With the external potential V(r) time-independent, we can write the inverse fermion propagator as the sum of a time-derivative and a position-dependent operator B(r), which we shall refer to as the Bogoliubov-de Gennes (bdg) operator: and h(r) is the one-particle Hamiltonian (2.12) Since B(r) is Hermitian, its eigenvalues are real and come in pairs.Moreover, the eigenvalues of the piece proportional to the identity in G −1 0 are given by iℏω n , where ω n (n ∈ Z) are the Matsubara frequencies (see previous footnote).It is possible to perform the sum over them explicitly and, in the zero-temperature limit, the result takes the simple form lim β→∞ The factor 1/2 in front of the trace has a direct physical interpretation, since we must distinguish between the "large" first-quantized Hilbert space, generated by all the modes (which is not a quantum system) and the "physical" Hilbert space of positive-energy modes.
One can obtain Eq. (2.13) by introducing the spectral zeta function associated with an operator : where the eigenvalues o α of  are labelled by a generic index α.With this, it is easy to show that For computational convenience, the last expressions above are cast in terms of the positive-definite operator B 2 .In order to compute such an object, it is useful to set up the heat kernel problem associated with B 2 .The heat kernel K B 2 (r 1 , r 2 ; τ) is a 2-by-2 matrix satisfying 5 which formally means that K B 2 (r 1 , r 2 ; τ) = ⟨r 1 |e −B 2 τ |r 2 ⟩.We denote the Dirac δ-function with a hat to avoid later confusion, and 1 2 is the 2-by-2 identity matrix.In turn, the coincident-point limit of the heat kernel allows one to compute the zeta function associated with B 2 through its Mellin transform where Tr now just stands for the matrix trace.Effectively, this turns the computation of the one-loop determinant -namely, the Tr log(|G −1 0 |) term -into the quantum mechanical problem of a particle with Hamiltonian B 2 .We will see in the following section that this is a good starting point for evaluating the functional determinant perturbatively in the limit of large chemical potential (or, equivalently, large charge), even when the system is inhomogeneous due to the presence of an external potential.
Before proceeding, it is useful to reformulate Eqs.(2.17) in dimensionless quantities, which will generically be denoted with a bar.As mentioned, we will eventually solve this problem perturbatively, but we need not commit to any approximation at this stage.The goal at this point is rather to express the problem in an exact form suggestive of 5.As usual, one should not confuse the proper time τ with the physical time.
a perturbative resolution.Consider the case of a generic potential V(r) that vanishes at the origin and that (classically) confines the particles to a finite region of space.Correspondingly, let R cl be the length scale associated with the smallest distance from the origin such that µ(r) = 0, i.e., where the particle density vanishes classically.For concreteness (and since in the end we want to compute conformal dimensions), one can consider the harmonic potential 6 V(r) = mω 2 2 r 2 , in which case, This length scale allows us to work in terms of dimensionless coordinates defined as Moreover, upon defining B(r) ≡ 1 ℏµ B(r) and where everything is dimensionless.Note that, explicitly, we have with σ(r) ≡ σ 0 (r) ℏµ , h(r) ≡ h(r) ℏµ = − ℏ 2mµR 2 cl △ r + V(r) − 1 and V(r) ≡ 1 ℏµ V(r).In the particular case of the harmonic potential, we simply have V(r) = r2 .

The bulk expansion
We have seen that the computation of the effective action at leading order in N reduces to evaluating the trace of the absolute value (or the square) of the bdg operator (cf.Eq. (2.13)).In the absence of a confining potential, this is a standard calculation [19,20].In the presence of a confining potential, the problem is no longer translationally invariant and the particles are confined to a spherical region at the edge of which the particle density drops to zero.
In this section, we show how to perturbatively compute the free energy in the presence of a potential when the particle number is large, thereby providing the first explicit verification of the predictions of nonrelativistic large-charge eft.The underlying 6.Though we do not need to commit to this choice for now.
controlling parameter of this perturbative computation is the particle density, which naturally breaks down close to the edge of the particle cloud.However, eft considerations not only tell us that this had to be anticipated, but also that "exotic" contributions enter the large-charge expansion of the free energy.While the setup in this section is not suited to capture these contributions, we discuss in Section 4 how to perform the matching between the expansion in the bulk of the cloud with the solution close to the edge, and find further agreement with the eft predictions.

Wigner coordinates and Moyal product
The lack of translational invariance can be addressed through the introduction of socalled mixed, or Wigner, coordinates [29,30].The idea is to first write B(r) as a bilocal operator via where and then perform a Fourier transform.Introducing relative and center-of-mass positions a bilocal Fourier transform of a function A(r 1 , r2 ) can be defined as and the inverse transformation is We choose to put a bar on the momentum as well to indicate that it is dimensionless.
In some sense, this transform allows one to disentangle the "microscopic" dynamics of the system (associated with rij ) from the "macroscopic" properties (associated with Rij ) resulting from the external potential [30].For the moment, ε is an arbitrary real parameter which will be assigned physical meaning soon.From this definition of the Fourier transform, it is easy to show the following general statement [31,32].If C(r 1 , r2 ) is related to two other bilocal functions A and B by where ★ is the Moyal product, defined by with what we shall call the k-Poisson bracket: Of course, the k = 1 case is just the normal Poisson bracket. 7 This formalism was originally introduced to describe quantum mechanics in phase space, and is adapted to solve our heat kernel problem, which takes the form where k( R, p; τ) denotes the Fourier transform of K B2 (r 1 , r2 ; τ), as defined above, and with the dimensionless phase-space Hamiltonian given by Note that the gap σ is only a function of R because it encodes pointwise interactions.Since, at leading order in N, we can identify 1 ℏβ S 0 with the grand-canonical potential Ω, its zero-temperature limit can be expressed in terms of the solution k( R, p; τ) to the above heat kernel problem as on top of which one needs to impose the gap equation.The charge is then given by which can be inverted to find µ as a function of Q, and the (zero-temperature) free energy is In the specific case of a critical system in a harmonic trap V(r) = mω 2 2 r 2 (i.e., V( R) = R2 and ε = ω/(2µ), cf.Eq. (2.19)), the free energy is related to the conformal dimension of the lightest operator of charge Q as [33,34]

Large-charge expansion of the heat kernel
Within the large-N expansion, this general construction is thus far exact, although the fact that the Moyal product can be expanded as an asymptotic series in ε suggests a perturbative calculation of the heat kernel problem.It is convenient to canonically normalize the momentum in the phase-space Hamiltonian in Eq. (3.14) by fixing which is small if the (dominating) scale associated with the potential is smaller than the scale set by the charge density.This is the large-charge limit because, as we will show, Q ∼ ε −3 .Our goal is to compute Ω(µ) as an expansion in ε, up to quadratic order, and thus to compute the corresponding free energy to this order. 8  8.As a reminder, the expansion to this order covers the dynamics of the effective action (1.1) to 'NLO', in the language of [12].On general grounds we expect the linear piece in ε to vanish because Ω(µ) is the the saddle This amounts to a semiclassical analysis of the Hamiltonian, Eq. (3.14), which is consistent as long as the potential term is bigger than the kinetic term.In the bulk region, the gap profile σ(r) does not vary rapidly on the scale defined by its own vev: (bulk condition). (3.21) In terms of dimensionless variables, this turns into which is satisfied as long as σ( R) and its derivative with respect to R are both of order one.More precisely, this means that 1 − V( R) must be bigger than some new parameter δ, which we will evaluate below.The construction is analogous to the standard Wentzel-Kramers-Brillouin (wkb) approximation, in which the expansion in ℏ is valid in the bulk away from the turning points.For the moment, we concentrate on the bulk region, where 1 − V( R) > δ (in boundary layer theory this is called the "outer region" [35]).We will turn to the edge of the cloud (the "inner region"), which necessitates a separate treatment, in the next section.
Having clarified the interval of validity of our approximation, we can use standard phase-space quantum mechanics.With our choice of ε (3.19), the one-particle Hamiltonian in phase space takes the form h( R, p) = p2 + V( R) − 1. (3.23) The dependence on ε has been reabsorbed by the Fourier transform in the same manner that ℏ does not appear in the phase-space Hamiltonian in quantum mechanics.This expression shows explicitly how in Wigner coordinates the contribution of the external potential, which depends only on the center-of-mass coordinate R, is separated from the kinetic part depending only on the (dual) relative coordinate p.
In the large-charge regime ε ≪ 1, it is then natural to make the following Ansätze for the gap and the heat kernel:

.25)
value of the function Ω(σ), and Ω NLO (µ) is the value of its first variation.
In the following, we will need to solve the saddle-point condition δΩ/δσ = 0.It is convenient to reformulate it in terms of Σ 0 by the chain rule .26)so that the saddle-point condition reads δΩ/δΣ 0 = 0.All the information about the saddle (i.e., the values of the functions Σ j ) is contained in the variations with respect to Σ 0 .The same information is also contained in the variations with respect to any of the Σ j , but this is suppressed in the ε expansion because In the expansion of Ω, the function Σ j appears linearly at order ε j : and, by the equation above, ΩNLO (Σ 0 ) = Ω ′ LO (Σ 0 ), which vanishes at the saddle.It follows that, apart from the leading order, the saddle-point value ⟨Ω⟩ at order j does not depend on the saddle value ⟨Σ j ⟩.
The phase-space Hamiltonian is quadratic in the variables R and p, and it is easy to express the product b .29) where we have identified, order by order, and With this, the heat kernel equation can be understood hierarchically.Expanding the Moyal product in powers of ε with the k-Poisson bracket notation introduced in Eq. (3.9), and dropping the arguments to avoid cluttering the notation, the order-ε n heat kernel problem becomes with initial conditions (ic) lim τ→0 K n = δ 0n and lim τ→0 Kn = 0. Note that each line contains a finite number of contributions and, for a given n, the two equations are decoupled, since B0 = 0. We investigate the order ε 0 contributions in the next section, and we will elaborate further on this system of equations when we reach subleading contributions.

Leading order
At leading order (lo), the Moyal product is just a pointwise product.This means that, at this order, σ commutes with the momentum and may therefore be treated as effectively constant.In other words, for very large µ, the position-dependent µ(r) can be regarded as a slowly varying function, and the computation is formally the same as the one without potential.In particular, by dimensional analysis and locality, the gap Σ 0 ( R) must be proportional to the effective chemical potential µ(r), namely and we must thereby reproduce the standard mean-field result for the Bertsch parameter in the absence of an external potential [36] (see also, e.g., [19]).It should soon become clear that this arises naturally in the lo analysis.The lo contributions are contained in the n = 0 (and therefore j = k = 0) terms in Eq. (3.32).The heat kernel equations at this order are with ic K 0 ( R, p; 0) = 1 and K0 ( R, p; 0) = 0. Thus, K 0 = e −B 0 τ and K0 = 0, and Eq.(3.15) at lo becomes (3.35) The experimentally tunable s-wave scattering length mentioned at the beginning is simply given by a s = mu 2πℏ 2 , and it diverges at the Feshbach resonance, indicating that the system has reached its critical (or unitary, in this context) point [8].In our scheme it is not necessary to renormalize the coupling of Σ 2 0 , and criticality simply corresponds to the limit u 0 → 0.
We want to express the result in the form of an integral over space, which is by construction well-defined only for 1 − V( R) < δ.At leading order, however, the result is well-defined up to R cl and we do not need any extra terms, which would otherwise lead to an unphysical scheme dependence.
If the integral over τ is performed first, we obtain integrals of the form which are well-defined for n > 1/2 and −3 < m < 4n − 5, but can be analytically continued to other values in terms of Gaussian hypergeometric functions 2 F 1 : At leading order we find where The saddle equation at this order is where we have used the recursion equation and admits the solution with y 0 ≈ 1.1622 . . .(see also Figure 2).This analysis confirms the predicted form of Σ 0 ( R), Eq. (3.33).At large N and for a large chemical potential µ, there is an explicit realization of the Cooper mechanism (see Eq. (2.4) and below): the Stratonovich field σ(τ, r) acquires a non-trivial vev which spontaneously breaks the U(1) symmetry.The system is then in a superfluid phase, which is the underlying assumption of the large-charge eft, allowing an explicit computation of the Wilsonian coefficients.
At the saddle we obtain c 0 = 4π I 0,0 (y 0 ) ≈ 0.0841 . . ., which, using Eq.(1.5), allows us to extract the value of the Bertsch parameter: In the absence of an external potential, this value was obtained via mean-field theory in the celebrated work [36], and agrees with the large-N computations in [19], in which the nlo-in-N correction was also computed.Next, recall that in Eq. (3.38), R 3 cl (2mµ/h) 32 = ε −3 ≫ 1 is the large-charge regime.In fact, the leading dependence of the charge on the chemical potential (at criticality) is given by potential V(r) = mω 2 2 r 2 (i.e., V( R) = R2 ) with ω ≪ µ, we find Using Eq. (3.18), the leading dependence of the conformal dimension of the large-charge operator on the charge Q is therefore given by

Subleading corrections
By inspection of Eq. (3.32), it is easy to show that the components of the heat kernel are of the form where K n,0 ( R, p) = δ 0n and Kn,0 ( R, p) = 0.That is, these components are polynomials in τ of degree (at most) n + 1, multiplied by an overall factor of e −B 0 τ.With this, the trace in Eq. (3.15) can be decomposed explicitly order by order in ε, and the integrals over momentum and proper time can be performed and expressed in terms of the functions I m,n .At criticality, Once the above saddle equation is accounted for, this recipe gives the order-ε n bulk contribution to the zero-temperature limit of the grand-canonical potential at criticality (and one still needs to account for the edge contributions).

Next-to-leading order
At next order in ε, i.e., n = 1 in Eq. (3.32), we use the fact that {B 0 , K 0 } = 0, so that the heat kernel equations at this order can be put into the form with K1 ( R, p; 0) = 0, K1 ( R, p; 0) = 0, ( where the hat notation means Thus, the first correction to the leading-order heat kernel is and we have B 1 = Σ * 0 Σ 1 + Σ * 1 Σ 0 , so we obtain The gap equation is which admits the single solution Σ 1 = 0. We find that the first corrections in ε to the gap and to the grand potential both vanish, in agreement with the eft prediction: (3.57)
The gap equation is obtained by varying with respect to Σ 0 : and, using the relations between the I m,n , we find the profile with y 1 ≈ −0.004347 . . ., y 2 ≈ −0.1608 . . . .(3.63)This is the first non-trivial correction to the profile of the Stratonovich field.We stress once more that this expression is only valid up to a distance δ from the edge at R = 1.We will discuss in the next section what happens closer to the edge of the cloud, where we will need to merge our small-gradient expansion with the perturbation theory of the boundary dynamics, subject to appropriate matching conditions.As for the effective action, we find 9 S = 4π ℏ 3 βm 3/2 I 0,0 (y 0 ) ∫ dr (ℏµ − V(r)) 5/2 + 5 64 From this, we can read off the Wilsonian parameters: The value for c 2 satisfies the bound of [12] on the transverse response function: The transverse response itself must be negative for a stable condensate, requiring that c 2 > 0.
The fact that 2c 1 + 3c 2 > 0 also places us in the regime were, in the zero-temperature, low-momentum limit, in the absence of an external potential or boundary dynamics, 1 → 2 phonon splitting is forbidden.
It is important to stress that the form that we find for the gap and the Lagrangian density is precisely the one expected based on locality and dimensional arguments: as long as we only look at the bulk, the problem has two scales, µ and ∇V, so any physical quantity of dimension ∆(G), in the limit of large µ, must take the form (3.67) 9. Notably, we can express the action in terms of rational coefficients of I 0,0 (y 0 ).

The edge expansion
In the previous section we obtained an asymptotic expansion of the bulk contributions to the free energy in terms of the small parameter ε ∼ ω/µ.These contributions match precisely the predictions from the bulk eft, and allow us to compute the corresponding Wilsonian coefficients at leading order in 1/N.We have noted, however, that this expansion is only valid in the bulk, up to some distance δ away from the point where the particle density vanishes at leading order in ε.This is the typical setup studied in boundary-layer theory [35].The domain of a differential equation is split into two (or more) regions.In each region, one finds an asymptotic expansion of the solution, which may then be matched over an intermediate region where both approximations are valid.
In the case where exact solutions are known across multiple regions, one is instructed to patch these solutions smoothly at transition points.Crucially, this is to be distinguished from cases where only asymptotic expansions are known in different regions, wherein one must instead perform a matching.That is, given two such solutions supported in different regions, one needs to find an overlapping region (the intermediate limit) in which the two asymptotic expansions have the same functional form.
Being able to find such a matching imposes a constraint on the thickness δ of the boundary layer.Qualitatively speaking, the boundary layer has to be sufficiently thick that the corresponding asymptotic solution (controlled by δ) is rich enough to permit a matching with the functional form coming from the outer (bulk) solution.Concretely, δ is fixed by a dominant-balance argument in which two or more of the terms of the oneparticle Hamiltonian must be of the same order, allowing us to retain enough information to perform a matching.This identifies the so-called distinguished limit.For simplicity, here we will consider the case of a spherically-symmetric potential.Correspondingly, the boundary layer is a spherical shell around r ≡ |r| = 1, so we introduce the inner variable ū as r = 1 − ūδ.(4.1) In the limit of δ → 0, the metric is approximately flat in terms of a radial and two orthogonal directions.To see that, start from R 3 in polar coordinates, with z ∈ (−1, 1) and ϕ ∈ (0, 2π).Rescaling z = xδ and ϕ = ȳδ, so that x ∈ (−1/δ, 1/δ) and ȳ ∈ (0, 2π/δ), the metric becomes and the Laplacian is approximately The only possible dominant balance is obtained when the first two terms are of the same order, i.e., ε 2 δ 3 = (1). (4.6) The double-scaling limit ε → 0, δ → 0, such that ε 2 /δ 3 = 1, is our distinguished limit in which the operator h is where α = V′ (1) and, in particular, for the harmonic oscillator α = 2.This is the Airy operator, corresponding to the linearization of the confining potential at the turning point.
From the point of view of the gap profile, in contrast to the bulk condition in Eq. (3.21), the edge condition is that σ varies on the same scale as the one fixed by its own vev

.8)
In the edge region 10 , then, we expect to be able to write the Stratonovich field as an expansion in δ so that the condition becomes which is satisfied if σ1 edge ( ū) and its ū derivative are of order (1).The boundary condition for σedge is obtained by imposing that in a region where both the bulk and edge approximations are consistent, the known result from the bulk is seems to imply a continuous spectrum.However, we are ultimately computing the trace of the absolute value | B(1) | = ((B (1) ) 2 ) 1/2 , and (B (1) ) 2 is bounded from below and increases without bound in the u direction, so that its spectrum is manifestly discrete.Thus, we can use a standard zeta-function regulator and write the free energy as the Mellin transform of the heat kernel at coincident points where we have used the fact that the integral over b(1) is of order 1/δ 2 .In the case of the harmonic trap we find F edge ht = ℏ 4/3 m 1/3  µ 5/3 ω 2/3 + . . ., (4.20) which is precisely the expected scaling of the leading-order edge contribution from the eft [16].

Conclusions
In this work, we have studied the unitary Fermi gas in a trapping potential in 3 + 1 dimensions at zero temperature in the limit of large charge and large number of Fermion flavors N.This is a well-justified approximation at the microscopic level, providing a foundation for the Thomas-Fermi approximation at leading order.We present a clean algorithmic procedure for extracting data from perturbation theory, requiring no other inputs beyond dimensional and symmetry arguments and basic techniques of quantum mechanics.
Previous attempts have involved a variety of approximation techniques, which we put here on a sound theoretical basis: The large N limit turns the one-loop calculation into a justified approximation at N = ∞ (which is typically extrapolated to N = 1 in mean field theory); the large-charge limit provides us with a small parameter proportional to the ratio of the coupling ω of the harmonic potential and the chemical potential ϵ = ω/(2µ), which allows us to perform a Moyal (gradient) expansion of the functional determinant resulting from integrating out the original microscopic dof.In this way we have found a consistent microscopic derivation that confirms the eft prediction for the form of the bulk operators, their contribution to the energy, and the form of the gap equation.We have also derived the first three Wilsonian coefficients at leading order in N and nnlo in ω/µ.Though the bulk expansion breaks down near the droplet edge, we were able to estimate the parametric dependenceof the leading-order edge contributions to the energy starting from first principles, confirming the eft prediction of [16].We reserve the precise determination of the numerical coefficient for future work.
Further characterizing the Wilsonian coefficients in the boundary theory constitutes a clear area for further investigation.Beyond this, one can also imagine expanding the reach of our theoretical understanding of the broader system by relaxing the limits of the regime studied herein: For instance, by including 1/N corrections, or by leaving the unitary point along the lines of [37,38], which would provide another controlled setting for the analysis discussed in [39], etc.Perhaps most importantly, it would be important to connect with experimental observations.As a first step, one might consider the doubly-integrated axial density [11], which can be computed based on the results herein.On the other hand, it would be valuable to obtain experimental data (of the density profile, perhaps) with sufficient resolution to probe the concrete predictions set forth here, particularly in the regime beyond the reach of the Thomas-Fermi approximation.

Figure 1 -
Figure 1 -Charge density as a function of the radial position in the harmonic trap for ε = ω/µ = 1/20.Shown are the leading bulk contribution (red) and the bulk contribution including the first correction (blue).The bulk approximation is consistent outside the δ-layer (shaded), where the edge expansion is needed.

Figure 2 -
Figure 2 -Spontaneous symmetry breaking.The effective action evaluated at the saddle has a minimum for a non-zero value of the parameter y that controls the expectation value of the gap.Shown is −I 0,0 (y), which appears in the effective action.