# Universal ground-state properties of free fermions in a d-dimensional trap

David S. Dean, Pierre Le Doussal, Satya N. Majumdar, Grégory Schehr

► To cite this version:

David S. Dean, Pierre Le Doussal, Satya N. Majumdar, Grégory Schehr. Universal ground-state properties of free fermions in a d-dimensional trap. EPL - Europhysics Letters, 2015, 112 (6), pp.60001 (1-6). 10.1209/0295-5075/112/60001 hal-01250783

HAL Id: hal-01250783

https://hal.science/hal-01250783

Submitted on 5 Jan 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
Universal ground state properties of free fermions in a d-dimensional trap

David S. Dean,¹ Pierre Le Doussal,² Satya N. Majumdar,³ and Grégory Schehr³

¹Univ. Bordeaux and CNRS, Laboratoire Ondes et Matière d'Aquitaine (LOMA), UMR 5798, F-33400 Talence, France

²CNRS-Laboratoire de Physique Théorique de l'Ecole Normale Supérieure, 24 rue Lhomond, 75231 Paris Cedex, France

³Univ. Paris-Sud, CNRS, LPTMS, UMR 8626, Orsay F-91405, France

We study the ground state properties of N spinless free fermions in a d-dimensional confining potential. We find that any n-point correlation function has a simple determinantal structure that allows us to compute several properties exactly for large N. We show that the average density has a finite support with an edge, and near this edge the density exhibits a universal (valid for a wide class of potentials) scaling behavior for large N. The associated edge scaling function is computed exactly and it generalizes the 1d result known from random matrix theory. In addition, we calculate the kernel (that characterizes any n-point correlation function) for large N and show that, when appropriately scaled, it depends only on dimension d, but has otherwise universal scaling forms, both in the bulk as well as at the edges. The edge kernel, for higher d, generalizes the Airy kernel in one dimension, well known from random matrix theory.

PACS numbers: 05.30.Fk, 02.10.Yn, 02.50.-r, 05.40.-a

Over the past few years, experimental developments in the construction of optical traps and cooling protocols have allowed the study of systems of confined ultracold atoms [1, 2]. Non-interacting fermions or bosons at low temperatures are of particular interest as they exhibit purely quantum effects. For example, Bose-Einstein condensation [1, 2] has been observed experimentally in several cold atom systems. There are also nontrivial quantum effects in non-interacting fermionic atoms, arising purely from the Pauli exclusion principle. A well studied example is a system of N spinless fermions in a one-dimensional harmonic potential, $V(x) = \frac{1}{2}m\omega^2x^2$ [3–12]. At zero temperature, $T = 0$, the many-body ground state wavefunction $\Psi_0(x_1, \dots, x_N)$ can be easily computed from the Slater determinant yielding [4, 10, 12]

$$|\Psi_0(x_1, \dots, x_N)|^2 = \frac{1}{z_N} \prod_{i<j} (x_i - x_j)^2 e^{-\alpha^2 \sum_{i=1}^N x_i^2}, \quad (1)$$

where $\alpha = \sqrt{m\omega/\hbar}$ has the dimension of inverse length and $z_N$ is a normalization constant. The squared wavefunction $|\Psi_0(x_1, \dots, x_N)|^2$, characterizing the quantum fluctuations at $T = 0$, can then be interpreted as the joint distribution of the eigenvalues of an $N \times N$ random matrix belonging to the Gaussian Unitary Ensemble (GUE) [13, 14]. Consequently, several zero temperature properties in this one-dimensional fermionic system have been computed analytically [5–12] using established results from random matrix theory (RMT). For example, the average density of fermions is given, for large N, by the celebrated Wigner semi-circular law [13, 14]:

$$\rho_N(x) \approx \frac{\alpha}{\sqrt{N}} f_W \left( \frac{\alpha x}{\sqrt{N}} \right); \quad f_W(z) = \frac{1}{\pi} \sqrt{2 - z^2} \quad (2)$$

with sharp edges at $\pm\sqrt{2N}/\alpha$. These sharp edges get smeared for finite but large N over a width $w_N \sim N^{-1/6}$ and the density near the edge (say, the right one) is de-

<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) (a) Plot of $\rho_N(x, y)$ vs. $x$ and $y$ for $N = 28$ fermions in a harmonic potential in $d = 2$ obtained from the exact evaluation [15] of Eqs. (8) and (9). (b) Plot of $\rho_N(x, y)$ in the $(x, y)$ plane predicted from the asymptotic formula Eq. (4). (c) Comparison of (a) and (b) for $y = 0$ (the region below the asymptotic result is shown via the shaded region).</figcaption></figure>

scribed by a finite size scaling form [16, 17]

$$\rho_N(x) \approx \frac{1}{N w_N} F_1 \left[ \frac{x - \sqrt{2N}/\alpha}{w_N} \right] \quad (3)$$

with $w_N = N^{-1/6}/(\alpha\sqrt{2})$ and $F_1(z) = [Ai'(z)]^2 - z[Ai(z)]^2$, where $Ai(z)$ is the Airy function. Far left from the right edge, using $F_1(z) \sim \sqrt{|z|}/\pi$ as $z \to -\infty$, one can show that the scaling form (3) smoothly matches with the semi-circular density in the bulk (2). The edge scaling function $F_1(z)$ has been shown [9] to be universal, i.e., independent of the precise shape of the confining potential $V(x)$. In addition, all n-point correlation functions at zero temperature can be expressed as determinants constructed from a fundamental quantity called the kernel (see later for its precise definition). For large N, and away from the edge, the appropriately scaled ker-
nel converges [9] to the universal sine-kernel form known from RMT [13, 14]. In contrast, near the edges, it approaches [9] the universal Airy kernel also well known from RMT [17, 18]. Recently, for one-dimensional traps, several of the zero temperature results for the density as well as the kernel have been extended to finite temperature [19].

In many experimental setups, the optical traps are actually in higher dimensions $d > 1$. For $d > 1$, unfortunately there is no obvious relation between free fermions in a trap and RMT. Consequently, calculating analytically the zero temperature properties of spinless free fermions in a $d$-dimensional confining potential is a challenging problem. It turns out that while the bulk properties can be estimated rather accurately using the local density approximation (LDA, also known as the Thomas-Fermi approximation) [20, 21], the density near the edge, where the Thomas-Fermi approximation breaks down, is much harder to compute. Similarly, one would also like to analytically compute the edge scaling of the kernel for $d$-dimensional traps. Another natural question is whether the universality observed in the one-dimensional case persists in higher dimensions. In this Letter, we compute analytically the density and the kernel for large $N$ for a $d$-dimensional confining trap. We demonstrate that, for large $N$, the edge density and the edge kernel, suitably scaled, are described by universal scaling functions, independent of the details of the trap potential.

It is useful to summarize our main results. For simplicity, we first focus on the harmonic potential $V(r) = \frac{1}{2}m\omega^2r^2$, for which we show that the average density of the fermions at zero temperature for large $N$ is given by

$$ \rho_N(\mathbf{x}) \approx \frac{1}{N} \left( \frac{m}{2\pi\hbar^2} \right)^{d/2} \frac{[\mu - \frac{1}{2}m\omega^2r^2]^{d/2}}{\Gamma(d/2 + 1)}, \quad (4) $$

where $r = |\mathbf{x}|$ and $\mu \approx \hbar\omega[\Gamma(d+1)N]^{1/d}$ is the Fermi energy for large $N$. This is a generalization of the 1d semi-circular law found in RMT and is shown in a scaling form in Eq. (2) (see Fig. 1 for $d=2$). The density in Eq. (4) has a finite support of radius $r_{\text{edge}} = (\sqrt{2}/\alpha)[\Gamma(d+1)]^{1/(2d)}N^{1/(2d)}$ where $\alpha = \sqrt{m\omega/\hbar}$. For finite but large $N$, the density near this edge $\rho_N(\mathbf{x}) \equiv \rho_{\text{edge}}(\mathbf{x})$ gets smeared over a scale of order $N^{-1/(6d)}$ and is described by a finite size scaling form

$$ \rho_{\text{edge}}(\mathbf{x}) \approx \frac{1}{N} \frac{1}{w_N^d} F_d \left( \frac{r - r_{\text{edge}}}{w_N} \right), \quad (5) $$

where $w_N = b_d N^{-1/6d}$, with $b_d = [\Gamma(d+1)]^{-1/6d}/(\alpha\sqrt{2})$, denotes the width of the edge regime. We also find that the scaling function $F_d(z)$ is given explicitly by

$$ F_d(z) = \frac{1}{\Gamma(\frac{d}{2} + 1) 2^{\frac{4d}{3}} \pi^{\frac{d}{2}}} \int_0^\infty du \, u^{\frac{d}{2}} \text{Ai}(u + 2^{2/3} z), \quad (6) $$

a plot of which is shown in Fig. 2 for $d=1, 2, 3$. In $d=1$ it reduces to the RMT result $F_1(z) = [\text{Ai}'(z)]^2 - z[\text{Ai}(z)]^2$

<figure><img src="image_3.png" /><figcaption>FIG. 2. (Color online) Plot of the scaling functions $F_d(z)$ in Eq. (6) for $d=1, 2, 3$ (top to bottom) for the density near the edge. The oscillatory structure of the scaling function becomes less pronounced as the dimension $d$ increases.</figcaption></figure>

mentioned earlier [16, 17]. Furthermore, going beyond the one-point function, i.e., the density, we show that the $n$-point correlation function in the ground state can be written as an $n \times n$ determinant whose entries are given by a two-point function called the kernel. We compute the asymptotic properties of this kernel explicitly both in the bulk [see Eq. (17)] as well as at the edges [see Eq. (20)] and derive the associated scaling functions in all dimensions. In $d=1$, they reduce to the well known Sine and the Airy kernel respectively. We then show that these $d$-dependent scaling functions associated with the bulk and the edge kernel as well as the edge density in Eq. (6) are universal, i.e., they are valid for a generic class of spherically symmetric potentials [22].

We start with $N$ spinless free fermions in a $d$-dimensional potential $V(\mathbf{x})$. The single particle eigenfunctions $\psi_k(\mathbf{x})$ satisfy the Schrödinger equation, $\hat{H}\psi_k(\mathbf{x}) = \epsilon_k\psi_k(\mathbf{x})$, where $\hat{H} = -\hbar^2/(2m)\nabla^2 + V(\mathbf{x})$ is the Hamiltonian and the energy eigenvalues $\epsilon_k$ are labelled by $d$ quantum numbers denoted by $k$. At zero temperature, the ground state many-body wavefunction can be expressed as an $N \times N$ Slater determinant, $\Psi_0(\mathbf{x}_1, \cdots, \mathbf{x}_N) = (1/\sqrt{N!}) \det[\psi_k(\mathbf{x}_j)]$, constructed from the $N$ single particle wavefunctions with energy up to the Fermi level $\mu$ [23]. For a sufficiently confining potential, $\mu$ generically increases with increasing $N$ [9, 20]. For example, for a $d$-dimensional harmonic oscillator $V(\mathbf{x}) \equiv V(r) = m\omega^2r^2/2$, $\mu \approx \hbar\omega[\Gamma(d+1)N]^{1/d}$. Using $\det(A)\det(B) = \det(AB)$, the squared many-body wavefunction can then be expressed as a determinant

$$ |\Psi_0(\mathbf{x}_1, \cdots, \mathbf{x}_N)|^2 = \frac{1}{N!} \det_{1 \le i,j \le N} K_\mu(\mathbf{x}_i, \mathbf{x}_j) \quad (7) $$

where the kernel $K_\mu(\mathbf{x}, \mathbf{y})$ is given by

$$ K_{\mu}(\mathbf{x}, \mathbf{y}) = \sum_{\mathbf{k}} \theta(\mu - \epsilon_{\mathbf{k}}) \psi_{\mathbf{k}}^{*}(\mathbf{x}) \psi_{\mathbf{k}}(\mathbf{y}). \quad (8) $$
Here $\theta(x)$ is the Heaviside theta function. As in $d=1$ [see Eq. (1)], the squared wavefunction in (7) can be interpreted as the joint probability density of $N$ points in a $d$-dimensional space. By integrating out the $N-n$ coordinates of $|\Psi_0(\mathbf{x}_1, \dots, \mathbf{x}_N)|^2$ in Eq. (7), one can show (see Supp. Mat. [24]) that the $n$-point correlation function can be expressed as an $n \times n$ determinant whose entries are given precisely by the kernel $K_\mu(\mathbf{x}, \mathbf{y})$ in Eq. (8). In particular, for $n=1$, the density $\rho_N(\mathbf{x})$ is given by

$$ \rho_N(\mathbf{x}) = \frac{1}{N} K_\mu(\mathbf{x}, \mathbf{x}). \qquad (9) $$

Thus, the knowledge of the kernel provides a complete description of the statistical properties of the ground state.

To compute the kernel, we first establish a simple but very general relation between the kernel and the propagator of the single particle quantum problem. Taking derivative of Eq. (8) with respect to (w.r.t.) $\mu$ and performing a Laplace transform w.r.t. $\mu$ of the resulting relation, one finds

$$ \int_0^\infty \frac{\partial K_\mu(\mathbf{x}, \mathbf{y})}{\partial \mu} e^{-\mu \frac{t}{\hbar}} d\mu = \sum_{\mathbf{k}} e^{-t \frac{\mathbf{k} \cdot \mathbf{\hat{\mathbf{\mathbf{k}}}}}{\hbar}} \psi_{\mathbf{k}}^*(\mathbf{x}) \psi_{\mathbf{k}}(\mathbf{y}) . \quad (10) $$

The right hand side of (10) is simply the single particle propagator $G(\mathbf{x}, \mathbf{y}; t) = \langle \mathbf{x} | e^{-\frac{t}{\hbar} \hat{H}} | \mathbf{y} \rangle$ in imaginary time. Integrating the left hand side of Eq. (10) by parts and inverting the Laplace transform, using Bromwich inversion formula, gives

$$ K_\mu(\mathbf{x}, \mathbf{y}) = \int_{\Gamma} \frac{dt}{2\pi i} \frac{1}{t} \exp\left(\frac{\mu}{\hbar}t\right) G(\mathbf{x}, \mathbf{y}; t), \quad (11) $$

where $\Gamma$ denotes the Bromwich integration contour. While Eq. (11) is general and holds for arbitrary potential $V(\mathbf{x})$, calculating the propagator $G(\mathbf{x}, \mathbf{y}; t)$ explicitly is hard for general $V(\mathbf{x})$. Below, we first focus on the specific case of a harmonic oscillator for which $G$ is known explicitly [25]

$$ G(\mathbf{x}, \mathbf{y}; t) = \left( \frac{\alpha^2}{2\pi \sinh(\omega t)} \right)^{d/2} e^{-\frac{\alpha^2}{2\sinh(\omega t)} Q(\mathbf{x}, \mathbf{y}; t)} \quad (12) $$

where $Q(\mathbf{x}, \mathbf{y}; t) = (\mathbf{x}^2 + \mathbf{y}^2) \cosh(\omega t) - 2 \mathbf{x} \cdot \mathbf{y}$. General potentials $V(\mathbf{x})$ will be considered later.

**Global density.** We first evaluate the global density $\rho_N(\mathbf{x})$ in Eq. (9) by putting $\mathbf{x} = \mathbf{y}$ in Eq. (11). The dominant contribution to the Bromwich integral in Eq. (11) with $\mathbf{x} = \mathbf{y}$ comes from the small $t$ region. Expanding the propagator to leading order for small $t$, the integral can be done explicitly to give the result in Eq. (4). Note that the normalization condition $\int d\mathbf{x} \rho_N(\mathbf{x}) = 1$ fixes the Fermi energy $\mu \approx \hbar\omega[\Gamma(d+1)N]^{1/d}$. The density thus has a radially symmetric finite support that vanishes at the edge as $\sim (r_{\text{edge}} - r)^{d/2}$ (see Fig. 1), where $r_{\text{edge}} = \sqrt{2\mu/(m\omega^2)} \sim N^{1/(2d)}$. Since $N$ particles are packed within a volume of radius $r_{\text{edge}} \sim N^{1/(2d)}$, the

typical inter-particle distance $\ell_{\text{typ}}$ can be estimated very simply: $N\ell_{\text{typ}}^d \sim r_{\text{edge}}^d$, implying $\ell_{\text{typ}} \sim N^{-1/(2d)}$.

**Edge density.** We next investigate the density near $r_{\text{edge}}$ for finite but large $N$. To derive the asymptotic edge behavior, we again start with the propagator $G(\mathbf{x}, \mathbf{y}; t)$, but now we set $r = |\mathbf{x}| = r_{\text{edge}} + z b_d N^{-\phi}$ where $\phi$ is yet to be determined and $b_d = [\Gamma(1+d)]^{-1/6d} / (\alpha\sqrt{2})$. Expanding the propagator for small $t$ and keeping terms up to order $O(t^3)$, we find that $\phi = 1/(6d)$ in order that the two leading terms scale in the same way for large $N$ with $z$ fixed (see Supp. Mat. [24]). Subsequently, evaluating the kernel in Eq. (11) and using Eq. (9), upon identifying $w_N = b_d N^{-1/(6d)}$, the edge density satisfies the scaling form in Eq. (5) where the scaling function is given by

$$ F_d(z) = (4\pi)^{-d/2} \int_{\Gamma} \frac{d\tau}{2\pi i} \frac{1}{\tau^{d/2+1}} e^{-\tau z + \tau^3/12} . \quad (13) $$

Using the integral representation of the Airy function, $Ai(z) = 1/(2\pi i) \int_{\Gamma} d\tau e^{-\tau z + \tau^3/3}$, the integral in Eq. (13) reduces to the expression announced in Eq. (6) (see Supp. Mat. [24]). This thus generalizes to arbitrary $d$ the $1d$ result, $F_1(z) = [Ai'(z)]^2 - z[Ai(z)]^2$, obtained from RMT [16, 17]. The asymptotic behaviors of $F_d(z)$ can be computed explicitly (see Supp. Mat. [24]) with the result

$$ F_d(z) \approx (8\pi)^{-\frac{d+1}{2}} z^{-\frac{d+3}{4}} e^{-\frac{4}{3}z^{3/2}} \quad \text{as } z \to \infty \quad (14) $$

$$ \approx \frac{(4\pi)^{-\frac{d}{2}}}{\Gamma(d/2 + 1)} |z|^{\frac{d}{2}} \quad \text{as} \quad z \to -\infty. \quad (15) $$

One can show that when $z \to -\infty$, i.e., when $r \ll r_{\text{edge}}$, the asymptotic behavior in Eq. (15) matches smoothly with the bulk density in Eq. (4).

**Bulk kernel.** We next consider the large $N$ scaling behavior of the kernel $K_\mu(\mathbf{x}, \mathbf{y})$ in Eq. (11) where the two points $\mathbf{x}$ and $\mathbf{y}$ are both far from the edge, while their relative separation $|\mathbf{x} - \mathbf{y}|$ is on the scale of the inter-particle distance $\ell_{\text{typ}} \sim N^{-1/(2d)}$. Once again, we start from our central equation (11). In the expression of the propagator $G(\mathbf{x}, \mathbf{y}; t)$ in Eq. (12), we first rewrite $Q(\mathbf{x}, \mathbf{y}; t) = (\mathbf{x} - \mathbf{y})^2 + (\mathbf{x}^2 + \mathbf{y}^2)(\cosh(\omega t) - 1)$. We then expand the propagator for small $t$ to leading order to obtain

$$ K_{\mu}(\mathbf{x}, \mathbf{y}) \approx \left( \frac{\alpha^2}{2\pi\omega} \right)^{\frac{d}{2}} \int_{\Gamma} \frac{dt}{2\pi i} \frac{1}{t^{\frac{d}{2}+1}} e^{\frac{(\mu - V(|\mathbf{x}|))t}{\hbar} - \frac{\alpha^2 (\mathbf{x} - \mathbf{y})^2}{2\omega t}} \quad (16) $$

where $V(|\mathbf{x}|) = V(r) = m\omega^2 r^2 / 2$. Fortunately, this integral can be done exactly (see Supp. Mat. [24]). We find that the bulk kernel has the scaling form, $K_\mu(\mathbf{x}, \mathbf{y}) \approx \ell^{-d} K_{\text{bulk}}(|\mathbf{x} - \mathbf{y}|/\ell)$, where $\ell = [N\rho_N(\mathbf{x})\gamma_d]^{-1/d}$ is the typical separation in the bulk and $\gamma_d = \pi^{d/2}[\Gamma(d/2 + 1)]$. The bulk scaling function is given explicitly by

$$ K_{\text{bulk}}(x) = \frac{J_{d/2}(2x)}{(\pi x)^{d/2}} \quad (17) $$

where $J_{d/2}(z)$ is the standard Bessel function of the first kind. In $d=1$, using $J_{1/2}(z) = \sqrt{2/(\pi z)} \sin z$, our result
in Eq. (17) again reduces to the standard sine-kernel in RMT [13]. The result in Eq. (17) is in full agreement with the heuristic derivation using the LDA [20] (see also Supp. Mat. [24]). However, the LDA becomes invalid near the edge where the local density is rapidly varying. We will see below that our approach yields exact results even in this edge regime where the LDA fails.

**Edge kernel.** Turning to the large $N$ behavior of the kernel $K_\mu(\mathbf{x}, \mathbf{y})$ near the edge, we set $\mathbf{x} = \mathbf{r}_{\text{edge}} + w_N \mathbf{a}$ and $\mathbf{y} = \mathbf{r}_{\text{edge}} + w_N \mathbf{b}$, following the scaling of the edge density in Eq. (5). Here $\mathbf{r}_{\text{edge}}$ denotes any point on the boundary of the support of the global density with $|\mathbf{r}_{\text{edge}}| = r_{\text{edge}} = \sqrt{2\mu/(m\omega^2)} \sim N^{1/(2d)}$. As before, the width $w_N = b_d N^{-1/(6d)}$ with $b_d = [\Gamma(d+1)]^{-1/(6d)} / (\alpha\sqrt{2})$. Thus $\mathbf{a}$ and $\mathbf{b}$ are dimensionless vectors. We substitute these scaling variables $\mathbf{x}$ and $\mathbf{y}$ in $Q(\mathbf{x}, \mathbf{y}; t) = (\mathbf{x}-\mathbf{y})^2 + (\mathbf{x}^2+\mathbf{y}^2)(\cosh(\omega t)-1)$ and expand $Q$ up to order $t^3$ for small $t$. Substituting these results in Eq. (11) and after a suitable change of variables (see Supp. Mat. [24]), one arrives at the edge behavior of the kernel

$$K_\mu(\mathbf{x}, \mathbf{y}) \approx \frac{1}{C_d w_N^d} \int_{\Gamma} \frac{d\tau}{2\pi i} \frac{1}{\tau^{\frac{d}{2}+1}} e^{-\frac{(\mathbf{a}-\mathbf{b})^2}{2^{8/3}\tau} - \frac{(a_n+b_n)\tau}{2^{1/3}} + \frac{\tau^3}{3}}, \quad (18)$$

with $C_d = (2^{4/3}\sqrt{\pi})^d$, and where $a_n = \mathbf{a} \cdot \mathbf{r}_{\text{edge}}/r_{\text{edge}}$ and $b_n = \mathbf{b} \cdot \mathbf{r}_{\text{edge}}/r_{\text{edge}}$ are projections of $\mathbf{a}$ and $\mathbf{b}$ in the radial direction. One can make a further simplification of Eq. (18) by using the integral representation of the diffusive propagator

$$\frac{e^{-\frac{(\mathbf{a}-\mathbf{b})^2}{4 D \tau}}}{(4\pi D \tau)^{\frac{d}{2}}} = \int \frac{d^d q}{(2\pi)^d} e^{-D q^2 \tau - i \mathbf{q} \cdot (\mathbf{a}-\mathbf{b})}. \quad (19)$$

We choose $D = 2^{2/3}$ and use this in Eq. (18). Using subsequently the integral representation of the Airy function $Ai(z)$ mentioned earlier, we arrive at the scaling behavior of the edge kernel, $K_\mu(\mathbf{x}, \mathbf{y}) \approx \frac{1}{w_N^d} K_{\text{edge}}(\frac{\mathbf{x}-\mathbf{r}_{\text{edge}}}{w_N}, \frac{\mathbf{y}-\mathbf{r}_{\text{edge}}}{w_N})$ where

$$K_{\text{edge}}(\mathbf{a}, \mathbf{b}) = \int \frac{d^d q}{(2\pi)^d} e^{-i \mathbf{q} \cdot (\mathbf{a}-\mathbf{b})} A i_1 \left( 2^{\frac{2}{3}} q^2 + \frac{a_n + b_n}{2^{1/3}} \right), \quad (20)$$

where $Ai_1(z) = \int_z^\infty Ai(u) du$. Putting $\mathbf{a} = \mathbf{b}$ in Eq. (20), followed by an integration by parts, one can check that $K_{\text{edge}}(\mathbf{a}, \mathbf{a})$ reduces to $F_d(|\mathbf{a}|)$ in Eq. (6). Also one can verify, after a few steps of algebra (see Supp. Mat. [24]), that in $d=1$ Eq. (20) reduces to $K_{\text{edge}}(a, b) = K_{\text{Airy}}(a, b) = (Ai(a)Ai'(b) - Ai'(a)Ai(b))/(a-b)$ is the standard Airy kernel [17, 18].

**General potential.** Having obtained these explicit results for harmonic oscillator potential, one naturally wonders to what extent these results are universal, i.e., hold for more general potentials $V(\mathbf{x})$. To investigate this issue of universality, let us first note that for general $V(\mathbf{x})$, our central equation (11) still holds, though the Fermi

energy $\mu$ as well as the propagator $G(\mathbf{x}, \mathbf{y}; t)$ depend on $V(\mathbf{x})$. The dependence of the Fermi energy $\mu$ on $N$ can be easily estimated for large $N$ using semi-classical approximation [9, 20]. In contrast, $G(\mathbf{x}, \mathbf{y}; t)$ is hard to compute for general $V(\mathbf{x})$. However, as we have seen in the case of a harmonic oscillator, only the small $t$ behavior of $G$ is needed to evaluate Eq. (11) for large $N$. Fortunately, the small $t$ expansion of $G$ for general $V(\mathbf{x})$ can be carried out using perturbation theory [26] (see Supp. Mat. [24] for a derivation). Using this expansion to leading order in $t$ in Eq. (11), we find (see Supp. Mat. [24]) that the global density is given by

$$\rho_N(\mathbf{x}) \approx \frac{1}{N} \left( \frac{m}{2\pi\hbar^2} \right)^{d/2} \frac{[\mu - V(\mathbf{x})]^{d/2}}{\Gamma(d/2 + 1)}, \quad (21)$$

which generalizes the harmonic oscillator result in Eq. (4). Note that this result in Eq. (21) coincides with the one obtained using Thomas-Fermi approximation [20]. Evaluating Eq. (11) near the edge, as in the harmonic oscillator case, we find that the edge density is again given exactly by Eq. (5), where only the location of the edge $\mathbf{r}_{\text{edge}}$ and the width $w_N$ depend non-universally on $V(\mathbf{x})$, but the scaling function $F_d(z)$ is universal, i.e., is the same for all spherically symmetric potential $V(\mathbf{x}) = V(r)$ [22]. The bulk kernel is also given as in Eq. (17) for the harmonic oscillator, the only $V(\mathbf{x})$-dependence appears in the scale factor $\gamma_d$ used in defining the pair of dimensionless vectors $\mathbf{u}$ and $\mathbf{v}$. In a similar fashion, we find that the edge kernel, appropriately centered and scaled, is also universal (see Supp. Mat. [24]).

**Conclusion.** We have studied the ground state properties of $N$ spinless free fermions in a $d$-dimensional confining trap and shown that in any dimension, the $n$-point correlation functions in the ground state have a determinantal structure with universal scaling properties for large $N$. This determinantal structure for the free fermions problem is well known in $1d$, due to its direct connection with RMT. In this Letter, we have shown that this structure persists even in higher dimensions, allowing us to derive several new exact results. Our results recover the bulk properties predicted by the heuristic LDA. However, near the edge where this approximation fails, our method predicts new universal exact results in all dimensions $d$. Our results can be extended to finite temperature [26] in all dimensions. Finally, it would be interesting to see whether one can measure these universal edge scaling functions in experiments on cold atoms.

We thank C. Salomon for useful discussions. We acknowledge support from PSL grant ANR-10-IDEX-0001-02-PSL (PLD) ANR grant 2011-BS04-013-01 WALK-MAT and in part by the Indo-French Centre for the Promotion of Advanced Research under Project 4604-3 (SM and GS).
[1] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. **80**, 885 (2008).

[2] S. Giorgini, L. P. Pitaevski, and S. Stringari, Rev. Mod. Phys. **80**, 1215 (2008).

[3] F. Gleisberg, W. Wonneberger, U. Schlder, C. Zimmermann, Phys. Rev. A **63**, 602 (2000).

[4] T. L. Einstein, Ann. Henri Poincaré **4**, Suppl. 2, S811 (2003).

[5] P. Calabrese, M. Mintchev, and E. Vicari, Phys. Rev. Lett. **107**, 020601 (2011); J. Stat. Mech. P09028 (2011).

[6] E. Vicari, Phys. Rev. A **85**, 062104 (2012).

[7] M. Campostrini and E. Vicari, Phys. Rev. A **82**, 063636 (2010).

[8] A. Angelone, M. Campostrini, and E. Vicari, Phys. Rev. A **89**, 023635 (2014).

[9] V. Eisler, Phys. Rev. Lett. **111**, 080402 (2013).

[10] R. Marino, S. N. Majumdar, G. Schehr, P. Vivo, Phys. Rev. Lett. **112**, 254101 (2014).

[11] I. Pérez-Castillo, Phys. Rev. E **90**, 040102(R) (2014).

[12] P. Calabrese, P. Le Doussal, S. N. Majumdar, Phys. Rev. A **91**, 012303 (2015).

[13] M. L. Mehta, *Random Matrices* (Academic Press, Boston, 1991).

[14] P. J. Forrester, *Log-Gases and Random Matrices* (London Mathematical Society monographs, 2010).

[15] For the $d$-dimensional harmonic oscillator the single particle wave functions $\psi_k(x) \propto e^{-\alpha^2 r^2} \prod_{j=1}^d H_{kj}(\alpha x_j)$ where $\alpha = \sqrt{m\omega/\hbar}$, $k_j$'s are non-negative integers and $H_k$ is the Hermite polynomial of degree $k$. Filling the $N$ lowest energy states (with degeneracy counted) we can enumerate the kernel in Eq. (8) explicitly and thus obtain the density from Eq. (9).

[16] M. Bowick and E. Brézin, Phys. Lett. B **268**, 21 (1991).

[17] P. J. Forrester, Nucl. Phys. B **402**(3), 709 (1993).

[18] C. A. Tracy and H. Widom, Commun. Math. Phys. **159**, 151 (1994).

[19] D. S. Dean, P. Le Doussal, S. N. Majumdar, G. Schehr, Phys. Rev. Lett. **114**, 110402 (2015).

[20] Y. Castin, *Basic theory tools for degenerate Fermi gases*, in Proceedings of the International School of Physics Enrico Fermi, Vol. 164: *Ultra-cold Fermi Gases*, edited by M. Inguscio, W. Ketterle, and C. Salomon, Varenna Summer School Enrico Fermi (IOS Press, Amsterdam, 2006).

[21] D. A. Butts, D. S. Rokshar, Phys. Rev. A **55**, 4346 (1997).

[22] Actually, the results hold for a wider class of potentials which are not necessarily spherically symmetric – details will be published elsewhere [26].

[23] For finite $N$ note that the ground state can be degenerate with the Fermi level $\mu$ partially occupied. However, for large $N$, this degeneracy does not affect the asymptotic properties and one can assume that all states up to $\mu$ are fully populated.

[24] see Supplemental Material which cites Refs. [27], [28] and [29].

[25] R. P. Feynman, A. R. Hibbs, *Quantum Mechanics and Path Integrals*, (McGraw- Hill, New York, 1965).

[26] D. S. Dean, P. Le Doussal, S. N. Majumdar, G. Schehr, in preparation.

[27] I. S. Gradshteyn and I. M. Ryzhik, *Table of Integrals, Series and Products*, edited A. Jeffrey and D. Zwilinger (Academic Press, Elsevier, 2007), 7th ed.

[28] O. Vallée and M. Soares, *Airy Functions and Applications to Physics*, (Imperial College Press, London, 2004).

[29] N. Makri and W. H. Miller, Chem. Phys. Letters **151** 1 (1988).