# Noninteracting fermions at finite temperature in a d-dimensional trap: Universal correlations

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. Noninteracting fermions at finite temperature in a d-dimensional trap: Universal correlations. Physical Review A: Atomic, molecular, and optical physics [1990-2015], 2016, 94 (6), pp.063622 (1-41). 10.1103/PhysRevA.94.063622. [hal-01421984]

HAL Id: hal-01421984

https://hal.science/hal-01421984

Submitted on 20 Jan 2017

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>
# Noninteracting fermions at finite temperature in a $d$-dimensional trap: Universal correlations

David S. Dean

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

Pierre Le Doussal

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

Satya N. Majumdar and Gregory Schehr

LPTMS, CNRS, Université Paris-Sud, Université Paris-Saclay, 91405 Orsay, France

We study a system of $N$ noninteracting spinless fermions trapped in a confining potential, in arbitrary dimensions $d$ and arbitrary temperature $T$. The presence of the confining trap breaks the translational invariance and introduces an edge where the average density of fermions vanishes. Far from the edge, near the center of the trap (the so-called “bulk regime”), where the fermions do not feel the curvature of the trap, physical properties of the fermions have traditionally been understood using the local density (or Thomas-Fermi) approximation. However, these approximations drastically fail near the edge where the density vanishes and thermal and quantum fluctuations are thus enhanced. The main goal of this paper is to show that, even near the edge, novel universal properties emerge, independently of the details of the shape of the confining potential. We present a unified framework to investigate both the bulk and the edge properties of the fermions. We show that for large $N$, these fermions in a confining trap, in arbitrary dimensions and at finite temperature, form a determinantal point process. As a result, any $n$-point correlation function, including the average density profile, can be expressed as an $n \times n$ determinant whose entry is called the kernel, a central object for such processes. Near the edge, we derive the large-$N$ scaling form of the kernels, parametrized by $d$ and $T$. In $d=1$ and $T=0$, this reduces to the so-called Airy kernel, that appears in the Gaussian unitary ensemble (GUE) of random matrix theory. In $d=1$ and $T>0$ we show a remarkable connection between our kernel and the one appearing in the $(1+1)$-dimensional Kardar-Parisi-Zhang equation at finite time. Consequently, our result provides a finite-$T$ generalization of the Tracy-Widom distribution, that describes the fluctuations of the position of the rightmost fermion at $T=0$, or those of the largest single-fermion momentum. In $d>1$ and $T \ge 0$, while the connection to GUE no longer holds, the process is still determinantal whose analysis provides a new class of kernels, generalizing the $1d$ Airy kernel at $T=0$ obtained in random matrix theory. Some of our finite-temperature results should be testable in present-day cold-atom experiments, most notably our detailed predictions for the temperature dependence of the fluctuations near the edge.

# I. INTRODUCTION

## A. Overview

Over the past few years, experimental developments in trapping and cooling of dilute Bose and Fermi gases have led to spectacular progress in the study of many-body quantum systems [1,2]. Even in the absence of interactions, bosons and fermions display collective many-body effects emerging purely from the quantum statistics [3–5]. For noninteracting fermions, which we focus on here, the Pauli exclusion principle induces highly nontrivial spatial (and temporal) correlations between the particles. Remarkably, the recent development of Fermi quantum microscopes [6–8] provides a direct access to these spatial correlations, via a direct *in situ* imaging of the individual fermions, with a resolution comparable to the interparticle spacing. It is thus important to have a precise theoretical description of these correlations in such fermionic gases.

In many experimental setups, including the aforementioned Fermi quantum microscopes, the fermions are trapped by an external potential. This trapping potential breaks the translational invariance and generically induces an edge of the Fermi gas. Indeed, beyond a certain distance from the

center of the trap, the average density of fermions vanishes. Far from the edge, i.e., close to the center of the trap, the properties of the Fermi gas are well described by standard tools of many-body quantum physics, like the local density (or Thomas-Fermi) approximation (LDA) [4,9]. However, this approximation breaks down close to the edge where the density vanishes, and where the fluctuations (both quantum and thermal) are large [10]. The purpose of this paper is to develop a general framework, which encompasses the LDA (and actually puts it on firmer ground) and provides an analytical description of the edge properties of the Fermi gas in any spatial dimension $d$, and at finite temperature $T$.

This framework, whose main results were recently announced in two short papers [11,12], is based on random matrix theory (RMT) [13,14] and, more generally, on the theory of determinantal point processes [15–17]. The simplest example is the case of $N$ noninteracting spinless fermions in a one-dimensional harmonic potential $V(x) = \frac{1}{2}m\omega^2x^2$ at zero temperature $T=0$. In this case, by computing explicitly the ground-state wave function, one can show [13,18,19] that there exists a one-to-one mapping between the positions of the fermions $x_1, x_2, \dots, x_N$ and the (real) eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_N$ of random $N \times N$ Hermitian matrices with
independent Gaussian entries, the so-called Gaussian unitary ensemble (GUE) in RMT. Although this connection has certainly been known for a long time [13], it is only recently that the powerful tools of RMT have been used to compute the statistics of physical observables for fermions [18–22] (see Sec. III for an extended discussion of these applications). In particular, it is well known that, in the large-N limit, the (scaled) density of eigenvalues (or equivalently the density of fermions) has a finite support and is given by the celebrated Wigner semicircular law $f_W(z) = \frac{1}{\pi}\sqrt{2-z^2}$.

Going beyond the average density, the fluctuations at the edge of the Wigner sea have also generated a lot of interest in RMT. Of particular interest is the statistics of the largest eigenvalue $\lambda_{\max} = \max_{1 \le i \le N} \lambda_i$. Indeed, its probability distribution function (PDF), properly shifted and scaled (see Sec. III below for more detail) is given for large N by the celebrated Tracy-Widom (TW) distribution for GUE [23], which is now a cornerstone of extreme value statistics of strongly correlated variables. Since its discovery in RMT, this TW distribution (and its counterparts for other matrix ensembles) have emerged in a wide variety of systems, a priori unrelated to RMT. These include the longest increasing subsequence of random permutations [24], directed polymers [25,26], and growth models [27,28] in the Kardar-Parisi-Zhang (KPZ) universality class in (1+1) dimensions as well as for the continuum (1+1)-dimensional KPZ equation [29–32], sequence alignment problems [33], mesoscopic fluctuations in quantum dots [34], height fluctuations of nonintersecting Brownian motions over a fixed time interval [35,36], height fluctuations of nonintersecting interfaces in presence of a long-range interaction induced by a substrate [37], and also in finance [38] (see [39,40] for reviews). The TW distributions have been recently observed in experiments on nematic liquid crystals [41] and in experiments involving coupled fiber lasers [42]. From the aforementioned connection between fermions and RMT, it follows that the quantum fluctuations of the rightmost fermion $x_{\max}(T=0) = \max_{1 \le i \le N} x_i$ are also governed by the TW distribution for GUE. Hence, noninteracting spinless fermions in a one-dimensional harmonic trap at $T=0$ provide one of the simplest physical systems where the TW distribution appears.

It is natural to ask what happens at finite temperature $T > 0$ and/or in higher dimensions $d > 1$? These are natural and relevant questions both experimentally and theoretically. In this case, the connection with the Gaussian unitary ensemble of RMT is lost. Nevertheless, it is still possible to show that the system is a determinantal process (this is an exact statement at $T=0$ and any $d$ [12] and true only for large $N$ at $T > 0$ [11]). This means that all $n$-point correlation functions of the $d$-dimensional Fermi gas, $R_n(x_1, x_2, \dots, x_n)$, can be expressed as the determinant of an $n \times n$ matrix whose entries are given by $K_\mu(x_i, x_j)$ where the function $K_\mu(x, y)$ is called the kernel and depends on the chemical potential $\mu$. It is thus the central object of the theory as it gives access to the computation of all the correlation functions. Therefore, it is important to characterize the behavior of this kernel, in the limit of a large number of fermions $N \gg 1$, both in the bulk and at the edge of the Fermi gas. In two recent papers, we announced, giving few details, results for the limiting kernels first in the case $d=1$ and $T \ge 0$ in Ref. [11] and then in

$d > 1$, but only for $T=0$ [12]. In this paper, we present a detailed derivation of these results, which relies in particular on a path-integral representation of the kernel $K_\mu(x, y)$. This allows us to study the problem in any dimension $d \ge 1$ and for any finite temperature $T \ge 0$, hence generalizing our previous results [12] to finite temperature. We will show that this method allows us to recover, in a fully controlled way, the LDA results in the bulk but also allows one to compute the correlations at the edge. In addition, as we will demonstrate below, this path-integral method is extremely powerful as it allows us to treat a wide class of trapping potentials of the form $V(x) \sim |x|^p$ at large $|x|$, with $p > 0$, and demonstrate the universality of our results, both in the bulk and at the edge.

## B. Model

We study in this paper a model for $N$ spinless non-interacting fermions in a $d$-dimensional potential $V(x)$. It is described by an $N$-body Hamiltonian $\hat{H}_N = \sum_{j=1}^N \hat{H}_j$, where $\hat{H}_j = \hat{H}(p_j, x_j)$ is a one-body Hamiltonian of the form $\hat{H} = \hat{H}(p_j, x_j)$ with

$$ \hat{H} = \frac{\hat{p}^2}{2m} + V(\mathbf{x}), \quad \hat{p} = \frac{\hbar}{i} \nabla. \qquad (1) $$

We will consider here only confining potentials $V(x)$ such that the one-body Hamiltonian $\hat{H}$ admits an infinite number of bound states. One such class consists of confining potentials of the form $V(x)$ with $V(x) \sim |x|^p$ at large $|x|$, with $p > 0$, and here we will mainly focus here on this class (and only briefly mention some other cases). The simplest such confining potential is the isotropic $d$-dimensional harmonic oscillator

$$ V(x) \equiv V(r) = \frac{1}{2} m \omega^2 r^2, \quad r = |x|. \qquad (2) $$

## C. Outline and summary of the main results

The paper is organized as follows. In Sec. II, we provide the general framework to study the correlations of $N$ noninteracting fermions in an arbitrary confining potential at zero temperature $T=0$ and in arbitrary dimension $d$. The main result is the determinantal structure of the $n$-point correlations in the ground state, given in Eq. (15), with the associated kernel in Eq. (9). In Sec. III, we focus on the one-dimensional system ($d=1$), for the special case of a harmonic potential $V(x) = \frac{1}{2}m\omega^2x^2$, still restricted to $T=0$. This case is particularly interesting because of its connection, valid for any finite number of fermions $N$, with the Gaussian unitary ensemble (GUE) of random matrix theory (RMT) [see Eqs. (20) and (21)]. From this connection, many interesting properties, which are summarized in that section, can be obtained for the fermion problem. In particular, this relation with RMT shows that the statistics of the position of the rightmost fermion, at zero temperature, is governed by the celebrated Tracy-Widom distribution for GUE [see Eqs. (47) and (48)]. In Sec. IV, we study the case of $N$ noninteracting fermions in a one-dimensional ($d=1$) arbitrary potential $V(x)$ at finite temperature $T > 0$. The main result of this section is to show that the $n$-point correlations still have a determinantal structure in the limit of a large number of fermions $N \gg 1$ [see Eqs. (93) and (94)]. As explained in detail in that section,
this structure is due to a large extent to the equivalence, when $N \gg 1$, between the canonical and the grand-canonical ensemble. In Sec. V, we apply the general analysis performed in Sec. IV to the case of $N$ noninteracting fermions in a one-dimensional harmonic potential at $T > 0$. This case is of particular interest as it is exactly solvable and the limit of large $N$ can be studied in detail. The main results of this section concern the local correlations which are described, in the limit $N \to \infty$, by the limiting kernels given by Eq. (120) in the bulk (i.e., near the center of the trap) and by Eq. (134) at the edge of the Fermi gas. They generalize the well-known sine kernel [see Eq. (35)] (in the bulk) and the Airy kernel [see Eq. (37)] (at the edge), which both play a fundamental role in RMT. The resulting kernel at the edge (134) allows us to compute the fluctuations of the rightmost fermion at finite temperature as a Fredholm determinant (137), which generalizes the Tracy-Widom distribution (48) for finite temperatures. Quite remarkably, exactly the same Fredholm determinant has appeared in the context of stochastic growth models in the KPZ equation [see Eqs. (305) and (307)]. This establishes an unexpected connection between free fermions at finite temperature and the KPZ equation at finite time $t$ [see Eq. (309)]. In Sec. VI, we study the case of $N$ noninteracting fermions in a harmonic potential, in arbitrary dimension $d$ and at zero temperature. Our analysis is based on a path-integral representation of the correlation kernel, given in Eqs. (168) and (169), which can then be analyzed in a very elegant way in the large-$N$ limit. The most interesting results of that section are certainly the expression of the density profile [see Eqs. (201) and (203)] and the kernel [see Eqs. (227) and (228)] at the edge, the latter being a generalization of the Airy kernel (37) for any finite $d$. In Sec. VII, we provide the full analysis of the correlations for a general $d$-dimensional soft potential, of the form $V(\mathbf{x}) \sim |\mathbf{x}|^p$ (for large $|\mathbf{x}|$ and $p > 1$), at finite temperature $T > 0$. Using again a path-integral representation of the kernel, we show that the local correlations both in the bulk and at the edge are universal, i.e., independent of the (smooth) confining potential considered here. The resulting universal correlation kernels in the bulk (274) and at the edge (296) generalize, respectively, the sine kernel (35) and the Airy kernel (37) for any finite dimension $d$ and temperature $T$. The last section (Sec. VIII) contains a discussion of our results, including the connection with the KPZ equation mentioned above, and our conclusions. Some technical details have been relegated in Appendixes A and B.

## II. CORRELATIONS FOR NONINTERACTING FERMIONS AT ZERO TEMPERATURE

### A. Many-body ground-state wave function

Let us start with $N$ fermions strictly at zero temperature. Consider first the single-particle eigenfunctions $\psi_k(\mathbf{x})$ which 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 labeled by $d$ quantum numbers denoted by $\mathbf{k} \in \mathbb{Z}^d$. Because of the confining potential, these quantum numbers labeled by $\mathbf{k}$ should not be identified with the usual momentum, which we denote by $\mathbf{p}$.

At zero temperature, the ground-state many-body wave function $\Psi_0$ can be expressed as an $N \times N$ Slater determinant

$$ \Psi_0(\mathbf{x}_1, \dots, \mathbf{x}_N) = \frac{1}{\sqrt{N!}} \det[\psi_{\mathbf{k}_i}(\mathbf{x}_j)]_{1 \le i, j \le N} \quad (3) $$

constructed from the $N$ single-particle wave functions labeled by a sequence $\{\mathbf{k}_i\}, i = 1, \dots, N$, with nondecreasing energies such that $\epsilon_{\mathbf{k}_i} \le \mu$ where $\mu$ is the Fermi energy. For a sufficiently confining potential, $\mu$ generically increases with increasing $N$ [4,18]. As an example, we consider the isotropic harmonic oscillator (2). In this case, the energy levels

$$ \epsilon_{\mathbf{k}} = \sum_{a=1}^{d} \left(k_a + \frac{1}{2}\right) \hbar \omega, \quad \psi_{\mathbf{k}}(\mathbf{x}) = \prod_{a=1}^{d} \phi_{k_a}(x_a), \\ \phi_k(x) = \left[ \frac{\alpha}{\sqrt{\pi} 2^k k!} \right]^{1/2} e^{-\frac{\alpha^2 x^2}{2}} H_k(\alpha x), \quad (4) $$

where the $k_a$'s are integers which range from 0 to $\infty$, $H_k(z)$ is the $k$th Hermite polynomial of degree $k$, and

$$ \alpha = \sqrt{m\omega/\hbar} \quad (5) $$

is the characteristic inverse length scale. Note that, in this example, each, non-ground-state, single-particle energy level is degenerate in $d > 1$. Hence, the $N$-body ground state is degenerate, whenever the last single-particle level is not fully occupied. This situation will be discussed later in Sec. IV. For now, since we are interested in the large-$N$ limit where this effect of degeneracy is subdominant, we will assume that the last level is fully occupied. In this case, for the harmonic oscillator, by filling up completely the levels up to $\mu$ one obtains $N = \sum_{\mathbf{k} \in \mathbb{Z}^d} \theta[\mu - \hbar\omega(k_1 + \dots + k_d)]$, where $\theta(x)$ is the Heaviside theta function. This leads for large $N$ to

$$ \mu \approx \hbar\omega[\Gamma(d+1)N]^{1/d}. \quad (6) $$

For more general potentials, the relation between $\mu$ and $N$ is usually more complicated (and detailed below), but given our assumption of an infinite number of bound states, we will always be able to study the limit of large $N \gg 1$, which is the subject of this paper.

### B. Quantum probability and determinantal structure of correlations

Consider now the quantum probability, i.e., the squared many-body wave function

$$ |\Psi_0(\mathbf{x}_1, \dots, \mathbf{x}_N)|^2 = \frac{1}{N!} \det[\psi_{\mathbf{k}_i}^*(\mathbf{x}_j)] \det[\psi_{\mathbf{k}_i}(\mathbf{x}_j)]. \quad (7) $$

Using the fact $\det(A^T)\det(B) = \det(AB)$, it can also be written as a determinant

$$ |\Psi_0(\mathbf{x}_1, \dots, \mathbf{x}_N)|^2 = \frac{1}{N!} \det K_\mu(\mathbf{x}_i, \mathbf{x}_j), \quad (8) $$

where we have defined the kernel $K_\mu(\mathbf{x}, \mathbf{y})$ as

$$ 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 (9) $$

As we will see later, this kernel will play a central role for the calculation of the correlations. For instance, one usually
defines the n-point correlation function $R_n(\mathbf{x}_1, \dots, \mathbf{x}_n)$ as

$$R_n(\mathbf{x}_1, \dots, \mathbf{x}_n) = \frac{N!}{(N-n)!} \int d\mathbf{x}_{n+1} \dots d\mathbf{x}_N |\Psi_0(\mathbf{x}_1, \dots, \mathbf{x}_N)|^2 \quad (10)$$

obtained by integrating over $N-n$ coordinates and keeping $n$ coordinates $\{\mathbf{x}_1, \dots, \mathbf{x}_n\}$ fixed. For $n=1$ this corresponds to the marginal density

$$R_1(\mathbf{x}) = N \int d\mathbf{x}_2 \dots d\mathbf{x}_N |\Psi_0(\mathbf{x}, \mathbf{x}_2, \dots, \mathbf{x}_N)|^2. \quad (11)$$

Incidentally, this is also related to the average local density of fermions

$$R_1(\mathbf{x}) = N\rho_N(\mathbf{x}), \quad \rho_N(\mathbf{x}) = \langle \hat{\rho}_N(\mathbf{x}) \rangle_0, \\ \hat{\rho}_N(\mathbf{x}) = \frac{1}{N} \sum_{i=1}^N \delta(\mathbf{x} - \mathbf{x}_i), \quad (12)$$

where $\langle \dots \rangle_0$ denotes the average w.r.t. the ground-state quantum probability in Eq. (7). The last equality follows from the indistinguishability of the fermions, i.e., the fact that the quantum probability $|\Psi_0(\mathbf{x}_1, \dots, \mathbf{x}_N)|^2$ is invariant under the exchange of any two coordinates. More generally, the $R_n$'s contain information about higher correlations of the local densities, e.g., one has

$$R_2(\mathbf{x}, \mathbf{y}) = \left\langle \sum_{1 \le i \ne j \le N} \delta(\mathbf{x} - \mathbf{x}_i) \delta(\mathbf{y} - \mathbf{x}_j) \right\rangle_0 \\ = N^2 \langle \hat{\rho}_N(\mathbf{x}) \hat{\rho}_N(\mathbf{y}) \rangle_0 - N \rho_N(\mathbf{x}) \delta(\mathbf{x} - \mathbf{y}) \quad (13)$$

and similarly for higher-order correlations.

Now we note that the kernel $K_\mu(\mathbf{x}, \mathbf{y})$ has the reproducing property

$$\int K_\mu(\mathbf{x}, \mathbf{z}) K_\mu(\mathbf{z}, \mathbf{y}) d\mathbf{z} = K_\mu(\mathbf{x}, \mathbf{y}), \quad (14)$$

which follows from the orthonormalization of the single-particle wave functions, $\int \psi_{\mathbf{k}_i}^*(\mathbf{z}) \psi_{\mathbf{k}_j}(\mathbf{z}) d\mathbf{z} = \delta_{\mathbf{k}_i, \mathbf{k}_j}$. If the kernel satisfies the reproducing property in Eq. (14), then there is a general theorem [13] that states that $R_n(\mathbf{x}_1, \dots, \mathbf{x}_n)$ in Eq. (10) can be expressed as an $n \times n$ determinant

$$R_n(\mathbf{x}_1, \dots, \mathbf{x}_n) = \det_{1 \le i, j \le n} K_\mu(\mathbf{x}_i, \mathbf{x}_j). \quad (15)$$

Note that this result (15) has been obtained here within the formalism of first quantization. It can also be derived within the formalism of second quantization: this is then a consequence of the Wick's theorem applied to fermionic (i.e., anticommuting) operators [3]. Indeed, the kernel in Eq. (9) can be expressed, in the second quantization formalism, as

$$K_\mu(x, y) = \langle \Psi_0 | \Psi^\dagger(x) \Psi(y) | \Psi_0 \rangle, \quad (16)$$

where $\Psi^\dagger(x)$ and $\Psi(y)$ are, respectively, the creation and the annihilation fermionic operators at positions $x$ and $y$ and $|\Psi_0\rangle$ is the ground state.

Any multiparticle probability distribution, whose n-point marginal can be expressed as the determinant of a kernel as in Eq. (15), will generally be referred to as a distribution with a *determinantal structure*. The associated random process corresponding to the random positions of the N fermions is then

called a d-dimensional *determinantal point process* [15,16]. Let us also point out a very simple consequence of this determinantal structure. Setting $n=1$ in Eq. (15) simply gives

$$R_1(\mathbf{x}) = K_\mu(\mathbf{x}, \mathbf{x}). \quad (17)$$

This implies, from Eq. (12), that the average density is given by the kernel evaluated at identical points

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

To summarize, the kernel $K_\mu(\mathbf{x}, \mathbf{y})$ is the key object for any determinantal process. Once we know the kernel, we can determine, in principle, any n-point correlation function by computing an $(n \times n)$ determinant [see Eq. (15)].

### III. ONE-DIMENSIONAL HARMONIC OSCILLATOR AT ZERO TEMPERATURE AND RMT

In this section, we consider the special case $d=1$ and the harmonic oscillator potential $V(x) = \frac{1}{2} m \omega^2 x^2$. In this special case, a host of analytical results for the zero-temperature quantum statistics have been derived over the years [18,43–47]. It turns out there is a close connection between the ground-state quantum probability $|\Psi_0(\mathbf{x}_1, \dots, \mathbf{x}_N)|^2$ of N fermions in a 1d harmonic trap and the joint probability distribution of N real eigenvalues of an $(N \times N)$ Gaussian Hermitian random matrix, known as the Gaussian unitary ensemble (GUE) in RMT. Although the connection between free fermions and GUE eigenvalues was known implicitly for a long time [13], this connection was first used explicitly, to our knowledge, in the context of studying the statistics of nonintersecting step edges on the vicinal surface of a crystal [48]. Subsequently, a connection between nonintersecting lines in presence of a potential $V(x) = x^2/2 + c/x^2$ (with $x \ge 0$ and $c > 0$) and the eigenvalue statistics of Wishart ensembles of RMT was established [37]. However, in the precise context of free fermions in a 1d harmonic trap, this connection to GUE eigenvalues was noticed and exploited only recently, first somewhat a posteriori in Ref. [18] and then more explicitly in Refs. [19,22], in the context of counting the number of fermions in an interval $[-L, L]$ in the ground state.

To establish the precise connection to GUE eigenvalues, we consider the first N single-particle levels with energies $\epsilon_k = (k+1/2)\hbar\omega$ where $k=0,1,2,\dots,N-1$. The many-body ground state is constructed by filling up these first N levels with N fermions. Thus, the ground-state energy is $E_0 = \frac{N^2}{2}\hbar\omega$, and the Fermi energy $\mu = (N-1/2)\hbar\omega$ corresponds to the highest occupied single-particle energy level in the many-body ground state. To construct the many-body ground-state wave function, we substitute the explicit single-particle harmonic oscillator wave functions [labeled by $k=0,1,2,\dots,(N-1)$ in Eq. (4)] in the $(N \times N)$ Slater determinant (3). This gives

$$\Psi_0(x_1, x_2, \dots, x_N) \propto e^{-\frac{\alpha^2}{2} \sum_{i=1}^N x_i^2} \det_{1 \le i, j \le N} [H_i(\alpha x_j)], \quad (19)$$

where we recall that $H_i(x)$ is the Hermite polynomial of degree $i$ and $\alpha = \sqrt{m\omega/\hbar}$ is an inverse length scale. By arranging the
rows and columns in the determinant, it is easy to see that it
can be reduced to a Vandermonde determinant up to an overall
constant $\det[H_i(\alpha x_j)] \propto \det[x_i^{j-1}] \propto \prod_{i<j}(x_i - x_j)$. Hence,
the ground-state quantum probability is given by [19]

$$
|\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 (20)
$$

where $z_N$ is a normalization constant.

Consider now a random (N × N) complex Hermitian matrix X with independent Gaussian entries, such that the joint distribution of independent matrix entries is given by Prob[X]dX ∝ exp[−Tr(X²)]dX. This joint distribution remains invariant under a unitary transformation X → U†XU, which justifies the name of such an ensemble of random matrices as the GUE [13]. Each realization of this matrix can be diagonalized to give N real eigenvalues {λ₁, λ₂, ..., λₙ} which are also random variables. What can be said about the joint distribution of these N eigenvalues? To obtain this eigenvalue distribution, one first makes a change of variables from the independent matrix entries to the eigenvalues and eigenvectors of X. Thanks to the rotational invariance, the eigenvector degrees of freedom decouple from the eigenvalues and hence can be integrated out. This provides an explicit expression for the joint distribution of eigenvalues [13]

$$
P(\lambda_1, \lambda_2, \dots, \lambda_N) = \frac{1}{Z_N} \prod_{i<j} (\lambda_i - \lambda_j)^2 e^{-\sum_{i=1}^N \lambda_i^2}, \quad (21)
$$

where $Z_N$ is the normalization constant. The Vandermonde
term $\prod_{i<j}(\lambda_i - \lambda_j)^2$ in Eq. (21) owes its origin to the
Jacobian of the change of variables from the matrix entries
to eigenvalues and eigenvectors [13].

Comparing Eqs. (20) and (21), it is clear that the quantum
statistics of the Fermion positions $x_i$'s in a 1d harmonic
trap at $T=0$ is identical, up to a trivial rescaling factor $\alpha$,
to the classical statistics of GUE eigenvalues $\lambda_i$'s [19]. The
routes of arrival to the identical joint distribution are, however,
quite different in the two problems. Interestingly, in the RMT
literature, to calculate further observables from this joint
eigenvalue distribution, the determinantal structure of this joint
distribution was noticed and exploited extensively [13,14].
In particular, the joint distribution in Eq. (21) was indeed

expressed as the determinant of a kernel,

$$
P(\lambda_1, \lambda_2, \dots, \lambda_N) = \frac{1}{N!} \det_{1 \le i,j \le N} K(\lambda_i, \lambda_j), \quad \text{with}
$$

$$
K(\lambda, \lambda') = \frac{1}{\sqrt{\pi}} e^{-\frac{1}{2}(\lambda^2 + \lambda'^2)} \sum_{k=0}^{N-1} \frac{1}{2^k k!} H_k(\lambda) H_k(\lambda'), \quad (22)
$$

where $H_k(z)$ is the $k$th Hermite polynomial. We recall that the Hermite polynomials are orthogonal on the real axis with respect to (w.r.t.) the Gaussian weight

$$
\int_{-\infty}^{\infty} e^{-z^2} H_n(z) H_m(z) dz = \sqrt{\pi} 2^n n! \delta_{n,m}. \quad (23)
$$

In the RMT literature, this orthogonal property of Her-
mite polynomials was exploited to derive the determinantal
structure of the joint distribution, hence is known in RMT
as the orthogonal polynomial method [13]. But, this is
precisely equivalent to the many-body quantum mechanics
of N fermions in a trap, once one recognizes that the joint dis-
tribution in Eq. (21) is just the square of the Slater determinant
and that $e^{-\lambda^2} H_k(\lambda)$ is just the $k$th single-particle eigenfunction
of a harmonic oscillator. The kernel $K(\lambda, \lambda')$ in Eq. (22) is also
identical, up to a trivial rescaling factor, to the one defined for
fermions in Eq. (9) in the case of a 1d harmonic oscillator.

Thus, to summarize, both the RMT and the many-body
free fermions essentially use the same method to analyze the
joint distribution. Hence, for this analysis, one does not really
need to know anything about random matrices. The starting
point is really the joint distribution in Eq. (20), which can then
be written as the determinant of a kernel $K_\mu(x,y)$. From the
general definition of the kernel in Eq. (9), upon summing up
the first N single-particle harmonic oscillator eigenfunctions
in Eq. (4), one obtains

$$
K_{\mu}(x, y) = \frac{\alpha}{\sqrt{\pi}} e^{-\frac{1}{2}\alpha^2(x^2 + y^2)} \sum_{k=0}^{N-1} \frac{1}{2^k k!} H_k(\alpha x) H_k(\alpha y). \quad (24)
$$

Note that because the Hermite polynomials are orthogonal
polynomials (23), they satisfy the Christoffel-Darboux iden-
tity [13,14] which allows us to perform the sum over *k* in
Eq. (24) explicitly, to yield (for *N* ≥ 1)

$$
K_{\mu}(x, y) = \begin{cases} \frac{e^{-\frac{\alpha^2}{2}(x^2+y^2)}}{\sqrt{\pi} 2^{N(N-1)!} \frac{H_N(\alpha x)H_{N-1}(\alpha y)-H_{N-1}(\alpha x)H_N(\alpha y)}{x-y}} & \text{for } x \neq y, \\ \frac{\alpha}{\sqrt{\pi}} \frac{e^{-\alpha^2 x^2}}{2^{N-1}(N-1)!} \{N[H_{N-1}(\alpha x)]^2 - (N-1)H_{N-2}(\alpha x)H_N(\alpha x)\} & \text{for } x = y. \end{cases} \tag{25}
$$

Consequently, the average density has the exact expression valid for any N:

$$
\rho_N(x) = \frac{1}{N} K_\mu(x, x) = \frac{\alpha}{N \sqrt{\pi}} e^{-\alpha^2 x^2} \sum_{k=0}^{N-1} \frac{1}{2^k k!} H_k^2(\alpha x) \quad (26)
$$

$$
= \frac{\alpha}{\sqrt{\pi}} \frac{e^{-\alpha^2 x^2}}{2^{N-1} N!} \{N[H_{N-1}(\alpha x)]^2 - (N-1)H_{N-2}(\alpha x)H_N(\alpha x)\}. \quad (27)
$$

A large number of precise analytical results for the ground-state quantum statistics of free fermions in a 1d harmonic trap have recently been predicted [18–22], essentially using the determinantal structure of the joint distribution and the explicit kernel in Eq. (24).
<figure><img src="image_2.png" /><figcaption>FIG. 1. The average density, for large $N$, has the Wigner's semicircular form $\rho_N(x) \approx \frac{\alpha^2}{\pi N} \sqrt{\frac{2N}{\alpha^2} - x^2}$, where $\alpha = \sqrt{m\omega/\hbar}$. The typical separation between particles in the bulk scales as $\sim N^{-1/2}$, where as near the edge it is much larger $\sim N^{-1/6}$.</figcaption></figure>

While the results in Eqs. (24)–(27) are exact for all $N$, it is useful and perhaps more interesting to see how the average density and the kernel behave asymptotically for large $N$. In the RMT literature, the large-$N$ asymptotics have been studied in great detail, mostly by using the asymptotic properties of the Hermite polynomials in Eqs. (24) and (27) [13,14]. For the benefit of readers not familiar with the RMT literature, we list below the principal RMT predictions for large $N$. For details of these derivations, the readers may consult Refs. [13,14].

### A. Large-$N$ RMT predictions for the average density

In the large-$N$ limit, the average density of fermions (equivalently that of GUE eigenvalues) is given by the celebrated Wigner semicircular 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 (28) $$

with sharp edges at $\pm\sqrt{2N}/\alpha$ (see Fig. 1). Note that the average density is normalized to unity, $\int \rho_N(x) dx = 1$ and $\rho_N(x)$ has the dimension of $\alpha$, i.e, inverse length. The result in Eq. (28) indicates that on an average there are more particles near the trap center $x = 0$ and less near the two edges. Thus, in the “bulk” of the Wigner sea, i.e., far away from the two edges, the density typically scales as $\rho_N(x) \sim \alpha N^{-1/2}$ for large $N$. This means that the typical interparticle separation in the bulk scales as $\ell(x) \sim 1/[N\rho_N(x)] \sim \frac{1}{\alpha} N^{-1/2}$ for large $N$.

In contrast, near the two edges, the particles are sparse (see Fig. 1) and the typical separation between two particles at the same edge scales as $\sim N^{-1/6}$ [13,14]. For finite but large $N$, the sharp edges at $\pm\sqrt{2N}/\alpha$ get smeared over a width $w_N \sim N^{-1/6}$. This is called the “edge” regime (see Fig. 1). The average density near the edge, for finite but large $N$, is described by a finite-size scaling form [49,50]

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

where we have set the width of the edge regime

$$ w_N = \frac{1}{\alpha \sqrt{2}} N^{-1/6}. \quad (30) $$

The scaling function is given exactly by [49,50]

$$ F_1(z) = [\mathrm{Ai}'(z)]^2 - z[\mathrm{Ai}(z)]^2, \quad (31) $$

where $\mathrm{Ai}(z)$ is the Airy function and $\mathrm{Ai}'(z)$ is its first derivative. The scaling function $F_1(z)$ has the asymptotic behavior

$$ F_1(z) \approx \begin{cases} \frac{1}{\pi} \sqrt{|z|} & \text{as } z \to -\infty, \\ \frac{1}{8\pi z} e^{-\frac{4}{3} z^{3/2}} & \text{as } z \to \infty. \end{cases} \quad (32) $$

Far to the left of the right edge, using $F_1(z) \sim \sqrt{|z|}/\pi$ as $z \to -\infty$ in Eq. (32), it is easy to show that the scaling form (29) smoothly matches with the semicircular density in the bulk (28). Recently, the edge scaling function $F_1(z)$ has been shown [18] to be universal, i.e., holds even for potentials different from the harmonic one, as long as it is smooth and confining.

### B. Large-$N$ RMT predictions for the kernel

The kernel in Eq. (24) can be analyzed similarly in the large-$N$ limit. For example, consider first the bulk with two points $x$ and $y$, both on the scale of the local interparticle separation $\ell(x)$, defined as

$$ \ell(x) = \frac{2}{\pi N \rho_N(x)}. \quad (33) $$

Taking the limit $N \to \infty$, and the separation between two points $|x - y| \to 0$, but keeping the ratio $|x - y|/\ell(x)$ fixed, it has been shown that the kernel $K_\mu(x, y)$ satisfies the scaling form

$$ K_\mu(x, y) \approx \frac{1}{\ell(x)} K^{\text{bulk}} \left( \frac{|x - y|}{\ell(x)} \right), \quad (34) $$

with the scaling function

$$ K^{\text{bulk}}(z) = \frac{\sin(2z)}{\pi z}. \quad (35) $$

This is the celebrated sine kernel which also turns out to be universal, i.e., independent of the precise shape of the trap potential $V(x)$ [18]. Note that when $z \to 0$, $K^{\text{bulk}}(z) \to 2/\pi$ and, consequently, the kernel $K_\mu(x, x) \to N \rho_N(x)$, in agreement with Eq. (18).

Similarly, near the edges (say the right edge at $x_{\text{edge}} = \sqrt{2N}/\alpha$), the kernel $K_\mu(x, y)$ in Eq. (24) can be similarly analyzed in the scaling limit: $N \to \infty$, $x \to x_{\text{edge}}$, $y \to x_{\text{edge}}$ but with the ratios $(x - x_{\text{edge}})/w_N$ and $(y - x_{\text{edge}})/w_N$ fixed. Here, $w_N$ denotes the width $w_N$ of the edge regime as defined in Eq. (30). In this scaling limit, one finds [13,14]

$$ K_\mu(x, y) \approx \frac{1}{w_N} K^{\text{edge}} \left( \frac{x - x_{\text{edge}}}{w_N}, \frac{y - y_{\text{edge}}}{w_N} \right), \quad (36) $$
where the two-variable scaling function is given by the so-called Airy kernel [13,14]

$$ 
\begin{aligned}
K^{\text{edge}}(a,b) &= K_{\text{Airy}}(a,b) = \frac{\text{Ai}(a)\text{Ai}'(b) - \text{Ai}'(a)\text{Ai}(b)}{a-b} \\
&= \int_{0}^{+\infty} du \, \text{Ai}(a+u)\text{Ai}(b+u). \quad (37)
\end{aligned}
$$

At coinciding points, it is easy to check that

$$ K_{\text{Airy}}(z,z) = F_1(z) = [\text{Ai}'(z)]^2 - z[\text{Ai}(z)]^2. \quad (38) $$

This fact, together with the definition $\rho_N(x) = K_\mu(x,x)/N$, then yields back the edge density result mentioned in Eq. (29).

### C. Statistics of the rightmost fermion, of fermion spacings, and of number fluctuations

From the connection with RMT, one can immediately obtain interesting predictions for various observables associated to 1d free fermions at $T = 0$. Indeed, in any determinantal point process [15,16], the full counting statistics can be obtained in terms of Fredholm determinants (denoted in this paper by Det). The Laplace transform of the probability $P_J(n)$ that there are exactly $N_J = n$ fermions in a given (arbitrary) subset $J$ of the real axis is given by

$$ \langle e^{-p N_J} \rangle = \text{Det}[I - (1 - e^{-p}) P_J K_\mu P_J], \quad (39) $$

where $P_J(x)$ is an indicator function, such that $P_J(x) = 1$ if $x \in J$ and $P_J(x) = 0$ otherwise, i.e., the projector on the interval $J$. In particular, the *hole probability*, i.e., the probability that there is exactly zero fermion in the subset $J$, is then

$$ P_J(n=0) = \text{Det}[I - P_J K_\mu P_J]. \quad (40) $$

There are important applications of this formula both in the bulk and at the edge, which we now discuss.

**Fermion spacing distribution.** One example in the bulk concerns the distribution of spacings (or gaps) between fermions, analog of the famous spacing distribution in RMT. Denoting by $\{x^{(i)}\}_{i=1,...,N}$ the ordered set of fermion positions (i.e., $x^{(1)} < x^{(2)} < ... < x^{(N)}$), we define the spacing $g = |x^{(i+1)} - x^{(i)}|$ between two consecutive fermions. The mean spacing near

the center (which we consider here) is $\bar{g} = 1/[N\rho_N(0)]$. The simplest guess for this spacing distribution is the famous Wigner surmise

$$ p_2^W(s) = \frac{32s^2}{\pi^2} e^{-4s^2/\pi}, \quad s = g/\bar{g} \quad (41) $$

which is normalized so that $\int_0^{+\infty} ds \, p_2^W(s) = \int_0^{+\infty} ds \, s \, p_2^W(s) = 1$, i.e., the mean fermion spacing is set to unity. As is well known, this is the exact result for $N = 2$. Thus, it is an exact statement for two fermions in a quadratic well (at $T = 0$). It also approximates rather well the *exact* distribution for a large number $N$ of fermions. The latter, close to the origin, can be obtained as

$$ p_2(g) = \frac{1}{N\rho_N(0)} \partial_g^2 D(g), \quad D(g) = \text{Det}[1 - P_{[0,g]} K_\mu P_{[0,g]}]. \quad (42) $$

In the bulk, setting $g = s/[N\rho_N(0)] = s^{\pi/2}\ell(0)$ (we made this rescaling by $\pi/2$ to conform to the standard convention used in RMT), one can replace the kernel in Eq. (42) by its limiting form, namely, the sine kernel in Eq. (35). The fermion spacing distribution is then described by the so-called Mehta-Gaudin distribution

$$ p_2(g) \approx N\rho_N(0)\tilde{D}''(s), \quad (43) $$

where $\tilde{D}(s)$ can be expressed in terms of a particular solution of a Painlevé V equation, denoted by $\sigma(x)$, such that

$$ \tilde{D}(s) = \exp \left[ \int_{0}^{\pi s} \frac{\sigma(x)}{x} dx \right], \quad (44) $$

where $\sigma(x)$ satisfies

$$ [x \sigma''(x)]^2 + 4[x \sigma'(x) - \sigma]\{x \sigma'(x) - \sigma(x) + [\sigma'(x)]^2\} = 0, \quad (45) $$

with the boundary condition $\sigma(x) \sim -x/\pi$ as $x \to 0$. From the Painlevé equation (45) the asymptotic behaviors of $\tilde{D}''(s)$ can be obtained as [13,51–53]

$$ \tilde{D}''(s) \sim \begin{cases} \frac{\pi^2 s^2}{3} - \frac{2\pi^4 s^4}{45} + \frac{\pi^6 s^6}{315} - \frac{\pi^8 s^8}{4050} - \frac{2\pi^{10} s^{10}}{14175} + \frac{11\pi^8 s^9}{496125} + \frac{2\pi^{10} s^{10}}{467775} + O(s^{11}), & s \to 0 \\ \frac{\pi^4}{16} \left(\frac{2}{\pi s}\right)^{1/4} \left[s^2 - \frac{2}{\pi^2} + o(1)\right] \exp\left[\frac{\ln 2}{12} + 3\zeta'(-1) - \frac{\pi^2 s^2}{8}\right], & s \to +\infty. \end{cases} \quad (46) $$

Note that from Eq. (43) one obtains that the average gap is given $\bar{g} = 1/[N\rho_N(0)] = \frac{\pi\ell(0)}{2}$, with $\ell(x)$ given in Eq. (33).

**Rightmost fermion statistics.** One important application of the formula (40) at the edge is as follows. In order to probe the statistics at the edge of the cloud of fermions, it is natural to consider the rightmost fermion $x_{\text{max}}(T=0) = \max_{1 \le i \le N} x_i$, where the quantum fluctuations of the positions of the $N$ fermions are described by the quantum probability in Eq. (20). Now, the cumulative probability distribution of $x_{\text{max}}(T=0)$ is precisely related to the hole probability in Eq. (40): Prob[$x_{\text{max}}(T=0) \le y$] is precisely the probability

that the interval $J = [y, +\infty)$ is free of particles. Using the expression (40) for the hole probability associated to the interval $J = [y = x_{\text{edge}} + s w_N, +\infty)$ (see below), one obtains that the typical quantum fluctuations of $x_{\text{max}}(T=0)$, correctly centered and scaled, are governed by the celebrated Tracy-Widom (TW) distribution for GUE, $F_2(x)$ [23]. Indeed, one has

$$ x_{\text{max}}(T=0) = x_{\text{edge}} + w_N \chi_2, \quad (47) $$

where, in the limit of large $N$, the cumulative distribution function (CDF) of the random variable $\chi_2$ is $F_2(s) = \text{Pr}(\chi_2 \le s)$,
which can be written as a Fredholm determinant [54]

$$ \mathcal{F}_2(s) = \mathrm{Det}(I - P_s K_{\mathrm{Airy}} P_s), \quad (48) $$

where $K_{\mathrm{Airy}}(a,b)$ is the Airy kernel given in Eq. (37) and $P_s$ is a projector on the interval $[s, +\infty)$. Note that $\mathcal{F}_2(s)$ can also be written in terms of a special solution $q(x)$ of the following Painlevé II equation [23]

$$ q''(x) = xq(x) + 2q^3(x), \quad q(x) \sim \mathrm{Ai}(x), \quad x \to \infty. \quad (49) $$

The TW distribution $\mathcal{F}_2(s)$ can then be expressed as

$$ \mathcal{F}_2(s) = \exp \left[ - \int_s^\infty (x-s) q^2(x) dx \right]. \quad (50) $$

In particular, its asymptotic behaviors are given by [55]

$$ \mathcal{F}_2(s) \sim \begin{cases} \tau_2 \frac{e^{-\frac{1}{2}|s|^3}}{|s|^{1/8}} [1 + \frac{3}{2^6|s|^3} + O(|s|^{-6})], & s \to -\infty \\ 1 - \frac{e^{-\frac{4}{3}x^{3/2}}}{16\pi s^{3/2}} [1 - \frac{35}{24s^{3/2}} + O(s^{-3})], & s \to +\infty \end{cases} \quad (51) $$

where $\tau_2 = 2^{1/24}e^{\zeta'(-1)}$ where $\zeta'(x)$ is the derivative of the Riemann zeta function.

**Number variance.** Another interesting observable is the number of fermions $N_L$ in a symmetric interval $[-L, L]$ (or in any fixed interval), which is also a random variable. Its mean is easily computed from the average density $\rho_N(x)$ in Eq. (12) as $\langle N_L \rangle = N \int_{-L}^{+L} \rho_N(x) dx$, which for large $N$ can be easily evaluated from the limiting semicircle law (28). What about the higher cumulants of this random variable, for instance, the variance $\mathrm{Var}(N_L)$, and eventually the full distribution of $N_L$? In RMT, the variance was computed a long time ago by Dyson [56], but only in the bulk limit when $L = s \ell(0)$ where $\ell(0)$ is the interparticle spacing in the center of the trap given in Eq. (33) for $x = 0$, and $s$ is a dimensionless number of order $O(1)$. In particular, for large $s$, one has

$$ \mathrm{Var}(N_L=s\ell(0)) \sim \frac{1}{\pi^2} \ln s + O(1), \quad s \to \infty \quad (52) $$

a result that can also be obtained using the LDA [4]. In the bulk regime one can show that the full distribution of $N_L$, properly centered and scaled, is a Gaussian [57,58]. It is only recently that the fluctuations of $N_L$, beyond the bulk regime, i.e., for $L \gg \ell(0)$, were studied. The variance $\mathrm{Var}(N_L)$ has indeed found a renewed interest [18,21,44–47,59] in the context of free fermions, thanks to its connection to entanglement entropy (see also below). Its dependence on $L$ was first studied numerically [18,45] for various trapping potentials, including the harmonic potential, and it was observed that it displays a striking nonmonotonic behavior. Recently, a full analytical computation of the variance as well as the full distribution of $N_L$ for fermions in a harmonic potential was performed in Refs. [19,22]. Using RMT tools, in particular Coulomb gas techniques, the variance was computed for any $L$ and large $N$, thus extending the analysis of Dyson [56] far beyond the bulk regime.

**Entanglement entropy.** The other interesting observable at $T=0$ is the Rényi entanglement entropy $S_q$ of the interval $[-L, L]$ around the trap center with the rest of the system. Consider the many-body fermionic system in its ground state,

so that the density matrix of the full system in this pure state is simply $\hat{\rho} = |\Psi_0\rangle\langle\Psi_0|$. Consider now the subsystem $A \equiv [-L, L]$ and let $\hat{\rho}_A = \mathrm{Tr}_{\tilde{A}}\hat{\rho}$ denote the reduced density matrix of the subsystem $A$, obtained by tracing out the complementary subsystem $\tilde{A}$ (so that $A$ and $\tilde{A}$ together constitute the full real line). The Rényi entanglement entropy of the subsystem $A$, parametrized by $q \ge 1$, is then defined as $S_q = \frac{1}{1-q} \ln \mathrm{Tr} \hat{\rho}_A^q$. In the limit $q \to 1$, this reduces to the standard von Neumann entropy $S_1 = -\mathrm{Tr}[\hat{\rho}_A \ln \hat{\rho}_A]$. For free fermions in a harmonic trap in one dimension at $T=0$, the Rényi entropy was studied numerically [18,45]. Exploiting the connection of the fermionic system to RMT, the Rényi entropy was recently computed analytically for large $N$, and for a wide range of $L$ [21]. For instance, it was shown [21] that for all $L$ such that $\sqrt{2N}/\alpha - L \gg w_N$ [where we recall that $w_N = 1/(\alpha\sqrt{2})N^{-1/6}$, see Eq. (30)], there is an exact relationship between the Rényi entropy and the number variance $\mathrm{Var}(N_L)$ discussed above:

$$ S_q = \frac{\pi^2}{6} \left(1 + \frac{1}{q}\right) \mathrm{Var}(N_L). \quad (53) $$

However, around the edge, this relation breaks down and computing the scaling behavior of $S_q$ near the edge remains a challenging open problem [21].

Let us close by indicating that, as a property of determinantal point processes which generalizes (39) (see, e.g., [15–17]), there exists a class of averages over the positions of the fermions which can be expressed exactly in terms of a Fredholm determinant (which is valid for any finite $N$)

$$ \left\langle \prod_{i=1}^{N} f(x_i) \right\rangle_0 = \mathrm{Det}(I - L_f), \\ L_f(x,y) = [1 - f(x)]K_\mu(x,y), \quad (54) $$

where $K_\mu(x,y)$ is the kernel associated with the determinantal process. In Eq. (54), the average is the quantum average in the ground state and the function $f(x)$ is arbitrary, provided the right-hand side exists. Such formula can be useful for studying linear statistics of free fermions [60].

### D. Momentum distribution function for fermions in a harmonic trap

To close this section, we point out the remarkable property that the quantum joint PDF in the ground state of $N$ noninteracting fermions in a harmonic potential is completely symmetric under the interchange of all positions with all momenta

$$ (\alpha x_1, \dots, \alpha x_N) \equiv_{\mathrm{in law}} \left( \frac{p_1}{\hbar \alpha}, \dots, \frac{p_N}{\hbar \alpha} \right) \\ \equiv_{\mathrm{in law}} (\lambda_1, \dots, \lambda_N)_{\mathrm{GUE}}, \quad (55) $$

each set being distributed as the eigenvalues of a GUE random matrix as indicated in formulas (20) and (21). This identity is true for arbitrary $N$. As a consequence, all properties obtained in this paper for the positions of the fermions are also valid for their momenta, as long as we are studying the harmonic confining potential. In particular, the density in momentum space $n(p)$ is also, at $T=0$, given by the Wigner
semicircle and exhibits a sharp edge at $p = p_{\text{edge}} = \hbar\alpha\sqrt{2N}$. The fluctuations of the largest momentum take the form, for large $N$,

$$p_{\text{max}} = \hbar\alpha\sqrt{2N} + \hbar\alpha^2 w_N \chi_2, \quad (56)$$

where $w_N$ is defined in Eq. (30) and $\chi_2$ is distributed according to the Tracy-Widom GUE distribution (50). All the predictions below for the correlations in the bulk and at the edge, at zero and finite temperature, and in any dimension, extend accordingly in momentum space.

## IV. CORRELATIONS FOR NONINTERACTING FERMIONS AT FINITE TEMPERATURE $T > 0$

We now want to study the effect of nonzero temperature for $N$ noninteracting fermions in an external confining potential. The discussion below holds for arbitrary dimensions, but we will focus for simplicity on the $d=1$ case. For the harmonic oscillator at $T=0$ one could use, as in the previous section, the techniques of RMT. However, at finite temperature, even for the harmonic oscillator potential, these direct RMT connections are lost and one needs to develop new techniques.

We will focus on the canonical ensemble at temperature $T = 1/\beta$ that corresponds to a fixed number of fermions $N$, which is often the situation studied in cold-atoms experiments. Before doing that, we first describe the energy basis of the Hilbert space, and the determinantal properties of the corresponding wave functions, in a more general setting. This general formalism is then applied to the harmonic oscillator, in the following section.

### A. $N$-fermion Hilbert space and occupation number basis

In order to study excited states, we need to consider the full Hilbert space of the $N$ particles. A natural basis of this Hilbert space is formed by the eigenstates of the $N$-particle Hamiltonian $\hat{H}_N$. For noninteracting fermions, these eigenstates and this basis can be constructed from the eigenstates of the single-particle Hamiltonian $\hat{H}$. In the case of the one-dimensional harmonic oscillator, the eigenfunctions $\phi_k(x)$ are given in Eq. (4). The associated energy eigenvalues are $\epsilon_k = \hbar\omega(k + 1/2)$ where $k$ is an integer which ranges from 0 to $\infty$. From these single-particle eigenstates, one can construct all many-body eigenfunctions of $\hat{H}_N$ by putting $N$ fermions in $N$ different single-particle levels indexed by $k_1 < k_2 < \dots < k_N$. The fermionic nature of the particles allows at most one particle in a given single-particle level. Hence, one introduces the set of occupation numbers, denoted by $\{n_k\}$, $k=0,1,2,\dots$ with $n_k = 0,1$, to label the many-body states, with $n_{k_1} = n_{k_2} = \dots = n_{k_N} = 1$ for the occupied single-particle states and $n_k = 0$ otherwise. They satisfy the constraint $\sum_{k=0}^{\infty} n_k = N$. The corresponding many-body eigenfunction is given by a Slater determinant, with the corresponding eigenenergy,

$$\Psi_{\{n_k\}}(\{x_i\}) = \frac{1}{\sqrt{N!}} \det_{1 \le i,j \le N} \phi_k(x_j), \quad E \equiv E_{\{n_k\}} = \sum_{k=0}^{\infty} n_k \epsilon_k \quad (57)$$

with, e.g., $E = \hbar\omega(k_1 + k_2 + \dots + k_N + N/2)$ for the harmonic oscillator.

An important property, already discussed and used in Sec. II for the ground state, but which extends to any $N$-body eigenstate, i.e., to all the excited states, is that the squared modulus of the wave function can be written as a determinant

$$\begin{aligned} |\Psi_{\{n_k\}}(x_1, \dots, x_N)|^2 &:= \frac{1}{N!} \left| \det_{1 \le i,j \le N} \phi_{k_i}(x_j) \right|^2 \\ &= \det_{1 \le i,j \le N} K(x_i, x_j; \{n_k\}), \end{aligned} \quad (58)$$

where the kernel $K(x,x'; \{n_k\})$ is indexed by the set of occupation numbers and is given by

$$K(x, x'; \{n_k\}) = \sum_{j=1}^{N} \phi_{k_j}^{*}(x) \phi_{k_j}(x') = \sum_{k=0}^{\infty} n_k \phi_{k}^{*}(x) \phi_{k}(x'). \quad (59)$$

Note that there is one kernel associated to each eigenstate of the energy operator $\hat{H}_N$. Furthermore, using the orthonormalization of the single-particle eigenfunctions one easily shows that these kernels obey the following property:

$$\int dz K(x,z; \{n_k\}) K(z,y; \{n'_k\}) = K(x,y; \{n_k n'_k\}) \quad (60)$$

for any two given sets $\{n_k\}$ and $\{n'_k\}$. Specializing (60) to $\{n_k\} = \{n'_k\}$ and using that $n_k^2 = n_k$ for any $k$, we see that each of these kernels satisfies the reproducing property (14). An immediate consequence, as in Sec. II, is that if the system is prepared in any of these states, the density and the correlations are given by determinants, as in Eqs. (18) and (15).

### B. Canonical measure and observables at finite $T$

Let us first recall the definition of the canonical partition function $Z_N(\beta)$ for $N$ noninteracting fermions at temperature $T = 1/\beta$ (in this paper we set the Boltzmann constant $k_B = 1$)

$$Z_N(\beta) = \text{Tr} e^{-\beta \hat{H}_N} = \sum_{k_1 < k_2 < \dots < k_N} e^{-\beta (\epsilon_{k_1} + \dots + \epsilon_{k_N})}, \quad (61)$$

where the sum is over all possible $N$-fermion eigenstates of $\hat{H}_N$ labeled as described above in terms of all possible distinct occupied single-particle eigenstates. It is also convenient to rewrite it using a labeling by occupation numbers as

$$Z_N(\beta) = \sum_{\{n_k\}} \left[ e^{-\beta \sum_{k \ge 0} n_k \epsilon_k} \delta\left(\sum_{k \ge 0} n_k, N\right) \right], \quad (62)$$

where $\sum_{\{n_k\}}$ denotes the sum over all the possible occupation numbers $n_k = 0,1$ for $k=0,1,2,\dots$. In Eq. (62), $\delta(i,j) = 1$ if $i=j$ and $\delta(i,j) = 0$ if $i \ne j$: this Kronecker delta function thus imposes the total number of particles to be exactly $N$, as we are working in the canonical ensemble.
In the canonical ensemble, the quantum joint probability distribution function of the positions $x_i$ of the fermions, $P_{\text{joint}}(x_1, \dots, x_N)$, is defined in terms of the N-body density matrix $\hat{\rho} = e^{-\beta \hat{H}_N}/Z_N(\beta)$ as

$$P_{\text{joint}}(x_1, \dots, x_N) = \langle x_1, \dots, x_N | \hat{\rho} | x_1, \dots, x_N \rangle = \frac{1}{Z_N(\beta)} \sum_{\{n_k\}} |\Psi_{\{n_k\}}(\{x_i\})|^2 e^{-\beta \sum_{k \ge 0} n_k \epsilon_k} \delta\left(\sum_{k \ge 0} n_k, N\right), \quad (63)$$

where the sum is over all many-body eigenstates. Using Eq. (57), we can rewrite the joint PDF of the particle positions in the canonical ensemble as the Boltzmann weighted sum of Slater determinants

$$P_{\text{joint}}(x_1, \dots, x_N) = \frac{1}{N! Z_N(\beta)} \sum_{\{n_k\}} \left| \det_{1 \le i,j \le N} \phi_{k_i}(x_j) \right|^2 e^{-\beta (\epsilon_{k_1} + \dots + \epsilon_{k_N})}, \quad (64)$$

where $Z_N(\beta)$ is the canonical partition function (61). It is easy to check that $Z_N(\beta)$ is such that the PDF $P_{\text{joint}}(x_1, \dots, x_N)$ is normalized to unity.

The first observable we want to compute is the mean density of fermions at finite temperature $T$ defined as

$$\rho_N(x) = \frac{1}{N} \sum_{i=1}^{N} \langle \delta(x - x_i) \rangle, \quad (65)$$

where from now on $\langle \dots \rangle$ means an average computed with the joint PDF (64). This amounts, up to a multiplicative constant, to integrating the joint PDF $P_{\text{joint}}(x, x_2, \dots, x_N)$ over the last $N-1$ variables. This amounts to the calculation of the following integral:

$$\rho_N(x) = \int_{-\infty}^{\infty} dx_1 \delta(x - x_1) \int_{-\infty}^{\infty} dx_2 \dots \int_{-\infty}^{\infty} dx_N P_{\text{joint}}(x_1, \dots, x_N), \quad (66)$$

where any of the two forms for $P_{\text{joint}}$ in Eqs. (63) and (64) can be inserted. More generally, we want to calculate the n-point correlation function $R_n(x_1, \dots, x_n)$ at temperature $T$ defined as

$$R_n(x_1, \dots, x_n) = \frac{N!}{(N-n)!} \int_{-\infty}^{\infty} dx_{n+1} \dots \int_{-\infty}^{\infty} dx_N P_{\text{joint}}(x_1, \dots, x_n, x_{n+1}, \dots, x_N). \quad (67)$$

The question is to handle these multiple integrals in the case of finite $T$. To understand the difficulty of this calculation, one can note that the joint PDF in Eq. (64) can be written as a determinant, as it is the case for $T=0$ [61]:

$$P_{\text{joint}}(x_1, \dots, x_N) = \frac{1}{N! Z_N(\beta)} \det_{1 \le i,j \le N} G(x_i, x_j, \beta\hbar) \quad (68)$$

in terms of the Euclidean propagator associated to the one-body Hamiltonian

$$G(x, y; t) = \langle y | e^{-\frac{i}{\hbar} \hat{H}} | x \rangle = \sum_{k} e^{-\frac{i}{\hbar} \epsilon_k} \phi_k^*(x) \phi_k(y). \quad (69)$$

Unfortunately, and at variance with the case of $T=0$, successive integrations over the coordinates $x_i$ do not preserve this determinantal structure. This is because the kernel inside the determinant no longer satisfies the reproducing property since

$$\int_{-\infty}^{\infty} dz \, G(x, z; \beta\hbar) G(z, y; \beta\hbar) = G(x, y; 2\beta\hbar) \quad (70)$$

which is clearly a different kernel. Hence, the evaluation of these integrals for arbitrary $N$ is very difficult.

### C. Saddle-point calculation and equivalence between canonical and grand-canonical ensembles

Fortunately, in the limit of large $N$, it is possible to use a saddle-point method to calculate the density and the n-point correlation functions (at fixed $n$). As we will see, this is a manifestation of the equivalence between the canonical and the grand-canonical ensembles for large $N$. Consider the density $\rho_N(x)$ in Eq. (66). Inserting there the expression (63) for $P_{\text{joint}}(x, x_2, \dots, x_N)$, and replacing $|\Psi_{\{n_k\}}(\{x_i\})|^2$ by the determinant of the kernel given in Eq. (58) we obtain

$$\rho_N(x) = \frac{1}{Z_N(\beta)} \sum_{\{n_k\}} \left[ \int_{-\infty}^{\infty} dx_1 \delta(x - x_1) \int_{-\infty}^{\infty} dx_2 \dots \int_{-\infty}^{\infty} dx_N \det_{1 \le i,j \le N} K(x_i, x_j; \{n_k\}) \right] e^{-\beta \sum_{k \ge 0} n_k \epsilon_k} \delta\left(\sum_{k \ge 0} n_k, N\right). \quad (71)$$
We now use the property of reproducibility of the kernel for each choice of {$n_k$} noted above in Eq. (60). From the theorem mentioned in Sec. II leading to Eqs. (15) and (18), we can rewrite this multiple integral as

$$ \rho_N(x) = \frac{1}{Z_N(\beta) N} \sum_{\{n_k\}} K(x, x; \{n_k\}) e^{-\beta \sum_{k \ge 0} n_k \epsilon_k} \delta \left( \sum_{k \ge 0} n_k, N \right), \quad (72) $$

where

$$ K(x, x; \{n_k\}) = \sum_{k=0}^{\infty} n_k |\phi_k(x)|^2. \quad (73) $$

Note that in the limit where $T \to 0$, the system is in the ground state characterized by $n_k = 1$ if $k = 0, 1, 2, \dots, N-1$ and $n_k = 0$ if $k \ge N$. Hence, in this limit, Eq. (72) reads as

$$ \rho_N(x) = \frac{1}{N} \sum_{k=0}^{N-1} |\phi_k(x)|^2. \quad (74) $$

To calculate the correlation functions given by the integral (67) we use the same method, and in particular the determinantal form (15) obtained after the $N-n$ integrations. We thus obtain the $n$-point correlations at finite temperature (67) as

$$ R_n(x_1, \dots, x_n) = \frac{1}{Z_N(\beta)} \sum_{\{n_k\}} \left[ \det_{1 \le i, j \le n} K(x_i, x_j; \{n_k\}) e^{-\beta \sum_{k \ge 0} n_k \epsilon_k} \delta \left( \sum_{k \ge 0} n_k, N \right) \right]. \quad (75) $$

To compute the expression in Eq. (75), we must evaluate the ratio of two sums over the occupation numbers {$n_k$}, each one constrained by $\sum_{k \ge 0} n_k = N$. Both in the numerator and the denominator, i.e., $Z_N(\beta)$ given by (61), we rewrite the constraint using an integral representation of the Kronecker delta symbol:

$$ \delta \left( \sum_{k \ge 0} n_k, N \right) = \int_0^{2\pi} \frac{d\lambda}{2\pi} \exp \left[ i\lambda \left( \sum_{k \ge 0} n_k - N \right) \right]. \quad (76) $$

Let us first consider the denominator $Z_N(\beta)$. It now reads as

$$ Z_N(\beta) = \int_0^{2\pi} \frac{d\lambda}{2\pi} e^{-iN\lambda} \sum_{\{n_k=0,1\}} e^{-\beta \sum_{k \ge 0} (n_k \epsilon_k - i\lambda n_k)}. \quad (77) $$

This representation thus makes the $n_k$ variables independent of each other and one can perform the sum over each $n_k$ separately. Note that each $n_k$ takes values 0 or 1. Performing the sum over $n_k$'s, for all $k$, we get

$$ Z_N(\beta) = \int_0^{2\pi} \frac{d\lambda}{2\pi} e^{-iN\lambda} \prod_k (1 + e^{-\beta(\epsilon_k - i\lambda)}). \quad (78) $$

Introducing the function

$$ J(\tilde{\mu}) = -\frac{1}{\beta} \sum_k \ln(1 + e^{-\beta \epsilon_k + \beta \tilde{\mu}}) \quad (79) $$

which is just the free energy in the grand-canonical ensemble at chemical potential $\tilde{\mu}$, we can rewrite the partition function simply as

$$ Z_N(\beta) = \int_0^{2\pi} \frac{d\lambda}{2\pi} e^{-\beta J(\frac{\tilde{\mu}}{\beta}) - i\lambda N}. \quad (80) $$

We now consider the full expression (75) and use (76) to rewrite it as

$$ R_n(x_1, \dots, x_n) = \frac{1}{Z_N(\beta)} \int_0^{2\pi} \frac{d\lambda}{2\pi} e^{-iN\lambda} \sum_{\{n_k\}} \det_{1 \le i, j \le n} \left[ \sum_{k=0}^{\infty} n_k \phi_k^*(x_i) \phi_k(x_j) \right] e^{-\beta \sum_{k \ge 0} (n_k \epsilon_k - i\lambda n_k)} \quad (81) $$

$$ = \frac{1}{Z_N(\beta)} \int_0^{2\pi} \frac{d\lambda}{2\pi} e^{-iN\lambda} \sum_{\{n_k\}} \sum_{k_1 < \dots < k_N} n_{k_1} \dots n_{k_N} e^{-\beta \sum_{k \ge 0} (n_k \epsilon_k - i\lambda n_k)} \left| \det_{1 \le i, j \le n} \phi_{k_i}(x_j) \right|^2. \quad (82) $$
In the second line we have used the (generalized) Cauchy-Binet-Andreief formula for determinants [62]:

$$ \det_{1 \le i, j \le n} \left[ \sum_{k=0}^{\infty} n_k \phi_k^*(x_i) \phi_k(x_j) \right] = \sum_{k_1 < \dots < k_N} n_{k_1} \dots n_{k_N} \left| \det_{1 \le i, j \le n} \phi_{k_i}(x_j) \right|^2 \quad (83) $$

valid for any set of {$n_k$}. To perform the sum over the occupation numbers, we introduce the notation for the “expectation value” of any observable $O[{n_j}]$ at fixed $\lambda$

$$ \langle O[{n_j}] \rangle_{\lambda} = e^{\beta J(\frac{i\lambda}{\beta})} \sum_{\{n_k\}} O[{n_j}] e^{-\beta \sum_{k \ge 0} (n_k \epsilon_k - i\lambda n_k)} \quad (84) $$

which is such that $\langle 1 \rangle_{\lambda} = 1$. This allows to rewrite (82) as

$$ R_n(x_1, \dots, x_n) = \frac{1}{Z_N(\beta)} \int_0^{2\pi} \frac{d\lambda}{2\pi} e^{-iN\lambda} e^{-\beta J(\frac{i\lambda}{\beta})} \sum_{k_1 < \dots < k_N} \langle n_{k_1} \dots n_{k_N} \rangle_{\lambda} \left| \det_{1 \le i, j \le n} \phi_{k_i}(x_j) \right|^2 \quad (85) $$

$$ = \frac{1}{Z_N(\beta)} \int_0^{2\pi} \frac{d\lambda}{2\pi} e^{-iN\lambda} e^{-\beta J(\frac{i\lambda}{\beta})} \sum_{k_1 < \dots < k_N} \langle n_{k_1} \rangle_{\lambda} \dots \langle n_{k_N} \rangle_{\lambda} \left| \det_{1 \le i, j \le n} \phi_{k_i}(x_j) \right|^2 \quad (86) $$

$$ = \frac{1}{Z_N(\beta)} \int_0^{2\pi} \frac{d\lambda}{2\pi} e^{-iN\lambda} e^{-\beta J(\frac{i\lambda}{\beta})} \det_{1 \le i, j \le n} \left[ \sum_{k=0}^{\infty} \langle n_k \rangle_{\lambda} \phi_k^*(x_i) \phi_k(x_j) \right] \quad (87) $$

and $\langle n_k \rangle_{\lambda} = \frac{1}{1+e^{\beta \epsilon_k - i\lambda}}$. In the second line we have used explicitly the independence of the variables $n_k$ at fixed $\lambda$, as seen from (84). In the third line we have used again the Cauchy-Binet-Andreief formula (in reverse). Thus, we finally obtain the n-point correlation function in the form

$$ R_n(x_1, \dots, x_n) = \frac{\int_0^{2\pi} \frac{d\lambda}{2\pi} [\det_{1 \le i, j \le n} K(x_i, x_j; \{n_k\}_{\lambda})] e^{-\beta J(\frac{i\lambda}{\beta}) - i\lambda N}}{\int_0^{2\pi} \frac{d\lambda}{2\pi} e^{-\beta J(\frac{i\lambda}{\beta}) - i\lambda N}} \quad (88) $$

with

$$ K(x, x'; \{n_k\}_{\lambda}) = \sum_{k=0}^{\infty} \langle n_k \rangle_{\lambda} \phi_k(x) \phi_k(x'), \quad \langle n_k \rangle_{\lambda} = \frac{1}{1 + e^{\beta \epsilon_k - i\lambda}} \quad (89) $$

and $J(\tilde{\mu})$ as defined in Eq. (79). At this stage, Eqs. (88), (89), and (79) provide an exact representation for the correlation function in the canonical ensemble for arbitrary $N$ and $n$, where the integrals over $\lambda$ still need to be performed. In particular, it holds for $R_1(x) = N\rho_N(x)$.

As a remark, we note that the crucial property which made this representation possible is the following identity:

$$ \left\langle \det_{1 \le i, j \le n} \left[ \sum_{k=0}^{\infty} n_k \phi_k^*(x_i) \phi_k(x_j) \right] \right\rangle = \det_{1 \le i, j \le n} \left[ \sum_{k=0}^{\infty} \langle n_k \rangle \phi_k^*(x_i) \phi_k(x_j) \right] \quad (90) $$

which holds for arbitrary averaging $\langle \dots \rangle$ for which the variables $n_k$ are independent. It can be proved using the Andreief formula twice, as done above. Here, we have further used the fact that the variables $n_k$ at fixed $\lambda$ are independent Bernoulli random variables. These identities have been used also in the mathematical literature [63] in the context of determinantal processes.

The next step is to evaluate the remaining integral over $\lambda$ using the saddle-point method, which becomes exact in the limit of large $N$. Let us first study the denominator in Eq. (88), i.e., the partition sum. There, the saddle point occurs, in our notation, at $\lambda = \lambda_{sp} = -i\beta\tilde{\mu}$ where $\tilde{\mu}$ is the chemical potential. It is related to the total number of particles $N$ as $N = -\partial J/\partial\tilde{\mu}$, which reads as

$$ N = \sum_{k=0}^{\infty} \frac{1}{e^{\beta(\epsilon_k - \tilde{\mu})} + 1} \quad (91) $$

Hence, the chemical potential $\tilde{\mu} = \tilde{\mu}(T, N)$ depends on $T$ and $N$. In the zero-temperature limit, as evident from the above equation, $\tilde{\mu}(T = 0, N) = \mu$, where $\mu$ is the Fermi level introduced in Sec. II.

In the thermodynamic language, this amounts to use the equivalence, in the large-$N$ limit, between the canonical ensemble and the grand-canonical ensemble. As is well known, the values of the average occupation numbers $n_k$ at the saddle point are given by the Fermi factor and denoted as

$$ \langle n_k \rangle := \langle n_k \rangle_{\lambda_{sp}}, \quad \langle n_k \rangle = \frac{1}{e^{\beta(\epsilon_k - \tilde{\mu})} + 1}, \quad N = \sum_{k=0}^{\infty} \langle n_k \rangle, \quad (92) $$

where we recall that for the harmonic oscillator $\epsilon_k = \hbar\omega(k + \frac{1}{2})$ for $k = 0, 1, 2, \dots$.
The same saddle-point analysis can be performed to calculate the correlation in Eq. (88). It remains valid as long as the quantity that is averaged does not grow too fast with $N$ [e.g., as $\exp(cN^a)$ where $a < 1$]. In this case, the value of the chemical potential at the saddle point is not modified. Here, this quantity is $\det_{1 \le i, j \le n} K(x_i, x_j; \{\langle n_k \rangle_\lambda\})$ and it is reasonable to expect that these conditions are satisfied. Therefore, in the large-$N$ limit, one obtains from Eq. (88) our main result

$$R_n(x_1, \dots, x_n) \approx \det_{1 \le i, j \le n} K_{\tilde{\mu}}(x_i, x_j), \quad (93)$$

where the finite-temperature kernel is given by

$$K_{\tilde{\mu}}(x, x') = \sum_{k=0}^{\infty} \frac{\phi_k^*(x)\phi_k(x')}{e^{\beta(\epsilon_k - \tilde{\mu})} + 1}, \quad (94)$$

and the chemical potential $\tilde{\mu}$ is fixed by Eq. (91). The case $n=1$ then yields the result for the density

$$\rho_N(x) \approx \frac{1}{N} \sum_{k=0}^{\infty} \langle n_k \rangle |\phi_k(x)|^2 = \sum_{k=0}^{\infty} \frac{|\phi_k(x)|^2}{e^{\beta(\epsilon_k - \tilde{\mu})} + 1}. \quad (95)$$

Note that we reserve the notations $\mu$ and $K_\mu(x, x')$ for the zero-temperature chemical potential and kernel, respectively, while we denote by $\tilde{\mu}$ and $K_{\tilde{\mu}}(x, x')$ their finite-temperature versions. When $T \to 0$, we recall that $\tilde{\mu} \to \mu$, but at finite $T > 0$, not only the chemical potential $\tilde{\mu}$ differs from $\mu$, but also the full kernel functions are different. Note also that in the limit $T \to 0$ the Fermi factor becomes $\theta(\mu - \epsilon_k)$ and the kernel becomes equal to the one associated with the ground state, given in Eq. (9) and the same result holds for the density.

Hence, we find that the correlations at fixed $n$ are asymptotically determinantal at large $N$ in the canonical ensemble. Alternatively, it is also possible to define the problem of fermions in an external potential directly in the grand-canonical ensemble. Indeed, the quantities $R_n(x_1, \dots, x_n)dx_1 \dots dx_n$ are the probabilities that there is a particle in each of the intervals $[x_i, x_i + dx_i]$, $1 \le i \le n$ (referred to as the correlation density in the mathematics literature). Clearly, such quantities exist and make sense for ensembles where the particle number varies. In fact in the grand-canonical ensemble, Eq. (93) is an exact equality. Hence, the kernel $K_{\tilde{\mu}}(x_i, x_j)$ exactly describes the statistics of a system in the grand-canonical ensemble at the chemical potential $\tilde{\mu}$ corresponding to $N$ for all values of $\tilde{\mu}$, and not only those corresponding to $N$ large. In the physics literature, this determinantal property of the grand-canonical ensemble of free fermions is usually derived using Wick's theorem [3,64] and has been known for a long time (see also Ref. [16]). Interestingly, in the mathematics literature, this property has been studied rigorously only rather recently, using Cauchy-Binet-Andreief identity in the context of general determinantal point processes [65]. Of course, both approaches provide identical results. In this paper, we have preferred the approach using Cauchy-Binet-Andreief identity to show that the determinantal property also holds in the canonical ensemble with fixed fermion number $N \gg 1$. We have shown that this is true provided the saddle-point solution exists. To prove the existence of such a saddle point rigorously is a challenging mathematical problem.

We end this section by commenting on the case where the single-particle spectrum is degenerate. In this case, the

many-body ground state may be degenerate, an example being the harmonic oscillator in $2d$. Each of these degenerate many-body ground states can be written as a Slater determinant and each of them constitutes a separate determinantal point process which can be studied along the lines of Sec. III. However, the zero-temperature limit of the finite-$T$ measure in Eq. (63) corresponds to taking a zero-temperature density matrix where each of these degenerate many-body ground states appears with equal probability. The resulting mixed state is not determinantal, as in the finite-$T$ case. However, using the equivalence between the canonical and the grand-canonical ensemble in the large-$N$ limit, this case can also be treated using Eq. (94) with $N$ given by (91), in the limit $\beta \to \infty$. Note that, everywhere in the formula given above, the sum over $k$ has to be understood as a sum over all possible single-particle states (including their degeneracies).

We now analyze these formulas (91)–(95) first in the bulk and then at the edge of the Fermi gas.

## V. HARMONIC OSCILLATOR IN ONE DIMENSION AT FINITE TEMPERATURE $T > 0$

Before applying the general formula for the finite-temperature kernel and correlations to the harmonic oscillator case, it is useful to discuss the relevant scales in the problem. We consider the harmonic oscillator potential $V(x) = \frac{1}{2}m\omega^2 x^2$ in $d=1$. We have seen that at $T=0$ there is a natural length scale associated to quantum fluctuations in the confining potential $1/\alpha = \sqrt{\hbar/m\omega}$. At $T=0$ there are two length scales $\ell(x) = \pi/N\rho_N(x)$ and $w_N = N^{-1/6}/(\alpha\sqrt{2})$ denoting, respectively, the interparticle distance in the bulk and the edge (see Fig. 1). A finite temperature introduces a length scale characterizing the width of a wave packet associated with a quantum particle, the thermal de Broglie wavelength $\lambda_T = \hbar\sqrt{2\pi/(mT)}$ obtained by the equating kinetic energy and $T$. Therefore, the thermal effects dominate over the quantum effects only if $\lambda_T$ is smaller than the typical interparticle distance, in which case the system behaves classically. In the bulk of the Fermi gas, comparing $\lambda_T$ and $l$ the quantum-to-classical crossover occurs at a temperature scale

$$T \sim N\hbar\omega. \quad (96)$$

Similarly at the edge, comparing $\lambda_T$ and $w_N$ we find that the corresponding crossover occurs at a much lower temperature

$$T \sim N^{1/3}\hbar\omega. \quad (97)$$

We will thus focus on these two regimes (96) and (97) in the following.

In addition, we know that the average density follows the Wigner semicircle law at zero temperature, a clear signature of the quantum effects. In the other limit of large temperature $T \gg N\hbar\omega$, the system behaves classically and we expect the standard Gibbs-Boltzmann distribution of independent particles

$$\rho_N(x) \xrightarrow{T \gg N\hbar\omega} \sqrt{\frac{\beta m\omega^2}{2\pi}} e^{-\frac{\beta}{2}m\omega^2 x^2}. \quad (98)$$

Our goal below is to study the quantum-to-classical crossover in the density as well as in the kernel.
A. High-temperature scaling and results in the bulk

As anticipated by Eq. (96), the bulk scaling regime corresponds to the limit $T \to \infty$, $N \to \infty$, while keeping fixed the following dimensionless variable:

$$y = N\hbar\omega/T = \beta\tilde{\mu}. \qquad (99)$$

Similarly, there is a length scale $\ell_T = \sqrt{2T/m\omega^2}$ associated with the high-temperature thermal fluctuations from Eq. (98), hence, we can consider the dimensionless length scale, also kept fixed

$$z = x/\ell_T = x\sqrt{m\omega^2/2T}. \qquad (100)$$

In this scaling limit, Eq. (91) fixing the chemical potential $\tilde{\mu}$ reads as

$$N = \sum_{k=0}^{\infty} \frac{1}{e^{\beta\hbar\omega(k+\frac{1}{2})}e^{-\beta\tilde{\mu}} + 1} \approx \int_{0}^{+\infty} dk \frac{1}{e^{ky/N}e^{-\beta\tilde{\mu}} + 1} = \frac{N}{y}\ln(1 + e^{\beta\tilde{\mu}}), \qquad (101)$$

where we could replace the sum by an integral since $\beta\hbar\omega \ll 1$. This yields the relation

$$e^{\beta\tilde{\mu}} = e^y - 1. \qquad (102)$$

Hence, in that scaling regime $\beta\tilde{\mu}$ is also of order $O(1)$.

1. Density of fermions in the bulk

We now analyze the density $\rho_N(x)$ given in Eqs. (95) and (92) which we evaluate for $x = z/\sqrt{\beta m\omega^2/2} = \sqrt{2N}/(\alpha\sqrt{y}) \gg 1$. After performing the change of variable $k = Np$ in the sum over $k$ in Eq. (95), one obtains

$$\rho_N(x) \approx \frac{1}{N} \sum_{k=0}^{\infty} \frac{\phi_k(x)^2}{e^{\beta\hbar\omega(k+\frac{1}{2})}e^{-\beta\tilde{\mu}} + 1} \approx \int_0^{\infty} dp \frac{\{\phi_{Np}[x = z\sqrt{2N}/(\alpha\sqrt{y})]\}^2}{e^{yp}(e^y - 1)^{-1} + 1}, \qquad (103)$$

where we have replaced $\tilde{\mu}$ by its value given in Eq. (102) and where $\phi_k(x)$ is given in Eq. (4). We now need an asymptotic expansion of $\phi_k(x)$ for large $k$ (and large $x$). This expansion is provided by the Plancherel-Rotach formula [as given for instance in Eqs. (3.10) and (3.11) of Ref. [66]]. For $-1 < X < 1$, one has

$$e^{-MX^2} H_M(\sqrt{2MX}) = \left(\frac{2}{\pi}\right)^{1/4} \frac{2^{M/2}}{(1-X^2)^{1/4}} M^{-1/4} (M!)^{1/2} g_M(X) \left[1 + O\left(\frac{1}{M}\right)\right] \qquad (104)$$

with

$$g_M(X) = \cos[MX\sqrt{1-X^2} + (M+1/2)\sin^{-1}X - M\pi/2]. \qquad (105)$$

Using these formulas (104) and (105) for $M = Np$ and $X = z/\sqrt{py}$ [see Eq. (103)] (taking into account that $X < 1$, i.e., $p > z^2/y$) one finds

$$\rho_N(x = z/\sqrt{\beta m\omega^2/2}) = \frac{\alpha\sqrt{2}}{\pi} \int_{z^2/y}^{\infty} dp \frac{1}{e^{yp}(e^y - 1)^{-1} + 1} (Np)^{-1/2} \frac{1}{\sqrt{1 - \frac{z^2}{py}}} [g_{Np}(z/\sqrt{py})]^2 \qquad (106)$$

$$= \frac{\alpha\sqrt{2}}{\pi\sqrt{y}\sqrt{N}} \int_{z^2}^{\infty} dq \frac{1}{e^{q}(e^y - 1)^{-1} + 1} \frac{1}{\sqrt{q - z^2}} [g_{Np}(z/\sqrt{q})]^2, \qquad (107)$$

where, in the second line, we have simply performed the change of variable $p = q/y$. To obtain the large-$N$ limit of Eq. (107) we notice that, thanks to the identity $\cos^2 x = 1/2 + \cos(2x)/2$, one can replace $[g_{Np}(z/\sqrt{q})]^2$, given in Eq. (105), in the integral over $q$ in Eq. (107) by $\frac{1}{2}$ (the remaining cosine being highly oscillating for large $N$ and thus subleading). If one finally performs the change of variable $q \to q + z^2$ in Eq. (107), one obtains finally

$$\rho_N(x = z/\sqrt{\beta m\omega^2/2}) = \frac{\alpha}{\pi\sqrt{y}\sqrt{2N}} \int_0^\infty dq \frac{1}{e^{q+z^2}(e^y-1)^{-1}+1} \frac{1}{\sqrt{q}} = -\frac{\alpha}{\sqrt{2N\pi y}} \text{Li}_{1/2}[-(e^y-1)e^{-z^2}], \qquad (108)$$

where $\text{Li}_n(x) = \sum_{k=1}^{\infty} x^k/k^n$ is the polylogarithm function. Hence, from Eq. (108) we obtain the main result of this section: the fermion density in the bulk takes the scaling form

$$\rho_N(x) \sim \frac{\alpha}{\sqrt{N}} R(y = \beta N \hbar \omega, z = x \sqrt{\beta m \omega^2 / 2}) \qquad (109)$$
with the bulk scaling function

$$R(y,z) = -\frac{1}{\sqrt{2\pi y}} \mathrm{Li}_{1/2}[-(e^y - 1)e^{-z^2}], \qquad (110)$$

which is plotted in Fig. 2. Note that it should not be confused with $R_n$ that denotes correlation functions.

We now show, from an asymptotic analysis of $R(y,z)$, that Eq. (109) interpolates between the Wigner semicircle (28) in the limit $T \to 0$ and the classical Gibbs-Boltzmann distribution for $T \to \infty$:

$$\rho_N(x) \underset{T \to \infty}{\longrightarrow} \sqrt{\frac{\beta m \omega^2}{2\pi}} \exp\left[-\frac{\beta}{2} m \omega^2 x^2\right], \qquad (111)$$

which holds also in the scaling limit $\beta \to 0$, $x \to \infty$ but keeping $x \sqrt{\beta}$ fixed (with the limit $N \to \infty$ already taken). Note that the physical mechanism behind this interpolation is very different from those found earlier in other matrix models [67,68].

To analyze the $T \to \infty$ and the $T \to 0$ limits of $\rho_N(x)$ in Eqs. (108) and (110), we need the following asymptotic behaviors of the polylogarithm function:

$$\mathrm{Li}_{1/2}(X) \sim X, \quad X \to 0 \qquad (112)$$

and

$$\mathrm{Li}_{1/2}(-e^X) \sim -\frac{2}{\sqrt{\pi}} X^{1/2}, \quad X \to \infty. \qquad (113)$$

From these behaviors in Eqs. (112) and (113), one finds the asymptotic behaviors of the scaling function $R(y,z)$ in Eq. (110):

$$R(y,z) \sim \begin{cases} \sqrt{\frac{y}{2\pi}} e^{-z^2}, & y \to 0 \\ \frac{\sqrt{2}}{\pi} \sqrt{1 - \frac{z^2}{y}}, & y \to \infty, z \to \infty \text{ with } z^2/y \text{ fixed.} \end{cases} \qquad (114)$$

From the first line of Eq. (114), one recovers the $T \to \infty$ limit where the density converges to the Gibbs-Boltzmann (Gaussian) form, given in Eq. (111). On the other hand, from the second line of Eq. (114), one obtains the $T \to 0$ limit of the density, which is given by the Wigner semicircle law (28).

## 2. Kernel and correlations in the bulk

We analyze the kernel $K_{\tilde{\mu}}(x,x')$ in the bulk where both $x = uN^{-1/2}/\alpha$ and $x' = u'N^{-1/2}/\alpha$ are close to the center of the trap (and $x - x'$ is of the order of the typical interparticle distance). Hence, we analyze the formulas (91) and (94) in the limit $N \to \infty$, $\beta \to 0$ keeping $y = \beta N \hbar \omega$ in Eq. (99) fixed. In this limit, the chemical potential $\tilde{\mu}$ is given by Eq. (102) and the large-$N$ analysis of $K_{\tilde{\mu}}(x,x')$ (94) can be performed along the same lines as done before for the density $\rho_N(x)$ yielding eventually Eq. (107). Indeed, using the Plancherel-Rotach asymptotic expansions (104) and (105), we

obtain

$$\begin{aligned} K_{\tilde{\mu}}(x = uN^{-1/2}/\alpha, x' = u'N^{-1/2}/\alpha) \\ \approx \frac{\alpha\sqrt{2N}}{\pi} \int_0^\infty \frac{dp}{\sqrt{p}} \frac{1}{e^y p (e^y - 1)^{-1} + 1} \\ \times g_{Np} \left( \frac{u}{N\sqrt{p}} \right) g_{Np} \left( \frac{u'}{N\sqrt{p}} \right). \end{aligned} \qquad (115)$$

From the explicit expression of $g_N(X)$ in Eq. (105), one obtains straightforwardly

$$g_{Np} \left( \frac{u}{N \sqrt{p}} \right) \approx \cos \left( 2u \sqrt{p} - Np \frac{\pi}{2} \right). \qquad (116)$$

Therefore, the product $g_{Np}(u/N\sqrt{p})g_{Np}(u'/N\sqrt{p})$ in Eq. (115) reads as, for large $N$,

$$\begin{aligned} g_{Np} \left( \frac{u}{N \sqrt{p}} \right) g_{Np} \left( \frac{u'}{N \sqrt{p}} \right) \\ \approx \frac{1}{2} \cos[2\sqrt{p}(u - u')] + \frac{1}{2} \cos[2\sqrt{p}(u + u') - Np\pi]. \end{aligned} \qquad (117)$$

The second term in Eq. (117) is highly oscillating in the large-$N$ limit and hence the leading contribution, once inserted in the integral over $p$ in Eq. (115), comes from the first term of Eq. (117), which is independent of $N$. Therefore, we finally

<figure><img src="image_3.png" /><figcaption>FIG. 2. Plot of the scaling function $R(y,z)$ associated with the density $\rho_N(x)$ (109), given in Eq. (110).</figcaption></figure>
obtain

$$
\begin{align}
K_{\tilde{\mu}}(x = uN^{-1/2}/\alpha, x' = u'N^{-1/2}/\alpha) & \nonumber \\
& \sim \alpha\sqrt{N} \frac{1}{\pi\sqrt{2}} \int_0^\infty \frac{dp}{\sqrt{p}} \frac{\cos[2\sqrt{p}(u-u')]}{e^{yp}(e^y-1)^{-1}+1}. \tag{118}
\end{align}
$$

Finally, performing the change of variable $p \rightarrow p/y$, one obtains the final form of the finite-temperature kernel in the bulk

$$
K_{\tilde{\mu}}(x, x') = \alpha N^{1/2} K_y^{\text{bulk}}[\alpha \sqrt{N}(x - x')], \quad (119)
$$

where

$$
K_y^{\text{bulk}}(v) = \frac{1}{\pi \sqrt{2y}} \int_0^{+\infty} dp \frac{\cos\left(\sqrt{\frac{2p}{y}}v\right)}{[1 + e^p/(e^y - 1)]\sqrt{p}} \quad (120)
$$

(see also Refs. [65,69] for alternative derivations of this
kernel). In the inset of Fig. 4 we show a plot of the two-
point correlation function in the bulk $g_y^{\text{bulk}}(v) = K_y^{\text{bulk}}(0)^2 -$
$[K_y^{\text{bulk}}(v)]^2$ for different scaled temperature parameter y.

B. Low-temperature scaling: Density and kernel at the edge

1. Density at the edge

We now focus on the density near the zero-temperature edge
at $x_{\text{edge}} = \frac{\sqrt{2N}}{\alpha}$. We recall that at zero temperature the density
strictly vanishes at the edge, and the edge region has a width
$w_N = \frac{1}{\alpha\sqrt{2}}N^{-\frac{1}{6}}$ which corresponds to the typical separation
between particles near the edge. To analyze how this density
profile gets modified near the edge we set

$$
x = \frac{\sqrt{2N}}{\alpha} + s w_N. \tag{121}
$$

We have also seen from Eq. (97) that in the edge region the crossover from quantum to classical regime occurs at temperature $T \sim N^{1/3}\hbar\omega$. Hence, we define the dimensionless (inverse temperature) parameter $b$,

$$
b = \frac{\hbar\omega}{T} N^{1/3}, \qquad (122)
$$

which will be kept fixed in the large-*N* and large-*T* limits
in this edge regime. From (99) we see that the variable *y* =
*b*N*²³/³ ≫ 1 in this regime. Hence, from (102) we can set βtilde*μ* ≃
*y* = *b*N*²³/³. We insert this value of βtilde*μ* in Eq. (95), and use ε*κ* =
*ħ*ω(*k* + 1/2). Making further a shift *k* − *N* = *m* (neglecting
the ½ factor compared to *N*) we obtain the following expression
for the density:

$$
\rho_N(x) \approx \frac{1}{N} \sum_{m=-N}^{\infty} \frac{[\phi_{N+m}(x)]^2}{\exp(bm/N^{1/3}) + 1}. \quad (123)
$$

Using the Plancherel-Rotach formula for Hermite polynomials
at the edge (see for instance Ref. [66]) yields

$$
\phi_{N+m} \left( \frac{\sqrt{2N}}{\alpha} + \frac{s}{\sqrt{2}\alpha} N^{-\frac{1}{6}} \right) \sim \sqrt{\alpha} \frac{2^{\frac{1}{4}}}{N^{\frac{1}{12}}} \text{Ai} \left( s - \frac{m}{N^{\frac{1}{3}}} \right), \quad (124)
$$

up to terms of order $O(N^{-2/3})$. Hence, by inserting this
asymptotic formula (124) into (123) and replacing the discrete

<figure><img src="image_4.png" /><figcaption>FIG. 3. Plot of the scaling function $F_{1,b}(s)$ for the density, given in Eq. (126), corresponding to two different (scaled) temperatures $b=0.5$ (dotted line) and $b=5$ (solid line). Inset: plot of $F_{1,b}(s)$ corresponding to $b \to \infty$ (zero temperature) shown here for comparison with the main plot.</figcaption></figure>

sum over $m \sim N^{1/3}$ by an integral, we obtain the scaling form
of the density near the edge

$$
\rho_N(x) \approx \frac{1}{N w_N} F_{1,b}\left(\frac{x - x_{\text{edge}}}{w_N}\right), \quad (125)
$$

where the finite-temperature scaling function $F_{1,b}(s)$ is ob-
tained as

$$
F_{1,b}(s) = \int_{-\infty}^{+\infty} du \frac{\mathrm{Ai}(s+u)^2}{1+e^{-bu}}. \quad (126)
$$

In the zero-temperature limit $b \to \infty$, the Fermi factor be-
comes a Heaviside step function, and we recover $F_{1,b=+\infty}(s) =
F_1(s)$ given in Eq. (31). In Fig. 3, we show how $F_{1,b}(s)$ behaves
for different values of the reduced inverse temperature $b$.
Note that the oscillations are more and more attenuated as
temperature increases.

*Density at the edge: Temperature dependence.* It is interesting to discuss the value of the density exactly at the edge $x = x_{\text{edge}}$. One has

$$
\begin{align}
N \rho_N(x_{\text{edge}}) &\approx N^{1/6} \alpha \sqrt{2} f_e(T/N^{1/3} \hbar \omega), \nonumber \\
f_e(t) &= t \int_{-\infty}^{+\infty} du \frac{\text{Ai}(ut)^2}{1+e^{-u}}. \tag{127}
\end{align}
$$

In the low-temperature scaling regime $T \sim N^{1/3}\hbar\omega$ the two
limiting behaviors are

$$
\rho_N(x_{\text{edge}}) \approx \begin{cases} \frac{\alpha \sqrt{2}}{3^{2/3} \Gamma(\frac{1}{3})^2} N^{-5/6}, & T \ll N^{1/3} \hbar \omega \\ \frac{(1-\sqrt{2}) \zeta(\frac{1}{2}) \alpha}{\sqrt{2\pi}} \frac{1}{N} \sqrt{\frac{T}{\hbar \omega}}, & T \gg N^{1/3} \hbar \omega. \end{cases} \quad (128)
$$

Note that we display the complete low-temperature series, as
an expansion in power of $T^2$, in Appendix B.

We recall that in the bulk high-temperature regime
$T \sim N\hbar\omega$, one has from (109) and (110)

$$
\rho_N(x_{\text{edge}}) \approx - \frac{\alpha}{\sqrt{2\pi N y}} \text{Li}_{1/2}(e^{-y} - 1), \quad y = N\hbar\omega/T \quad (129)
$$
which gives the two limiting behaviors

$$ \rho_N(x_{\text{edge}}) \approx \begin{cases} \frac{(1-\sqrt{2})\zeta(\frac{1}{2})\alpha}{\sqrt{2\pi}} \frac{1}{N} \sqrt{\frac{T}{\hbar\omega}}, & T \ll N\hbar\omega \\ \sqrt{\frac{m\omega^2}{2\pi T}}, & T \gg N\hbar\omega \end{cases} \quad (130) $$

which shows a perfect matching of the high-temperature end of the low-T (edge) scaling regime [second line of Eq. (128)], with the low-temperature end of the high-T (bulk) regime [first line of Eq. (130)]. We recall that in the high-T regime, the Fermi gas extends well beyond the T = 0 edge. Note that for fixed (large) N this density first increases as a function of T in the low-T regime and then exhibits a maximum for T ~ N in the high-T regime before decreasing again.

**Right tail of the density.** Let us consider the behavior of the density scaling function $F_{1,b}(s)$ to the right of the edge, for large positive $s$. The analysis of the integral in Eq. (126) in that limit was performed in Ref. [70]. It is found that there are two regimes depending on whether the parameter $\tilde{s} = s/b^2$ is smaller or larger than the critical value $\tilde{s}_c = \frac{1}{4}$:

$$ F_{1,b}(s) \approx \begin{cases} \frac{1}{4b^2\sqrt{s}\sin(2\pi\sqrt{s})} \exp\left(-\frac{4}{3}s^{3/2}\right), & 1 \ll s < \frac{b^2}{4} \\ \frac{1}{\sqrt{4\pi b}} \exp\left(-bs + \frac{b^3}{12}\right), & s > \frac{b^2}{4}. \end{cases} \quad (131) $$

Hence, we obtain a transition between a stretched exponential tail in the density, as in the zero-temperature case, to a pure exponential decay in the far tail for $s > b^2/4$. Thus, for a fixed value of the reduced temperature $b$ (not necessarily large), the decay is always exponential. Note that the preexponential factor exhibits a crossover from the two limiting cases indicated above, in the vicinity of $\tilde{s} = \frac{1}{4}$ [70].

**Left tail of the density.** For $s$ large and negative, the integral over $u$ in Eq. (126) is dominated by the region $u \in (-\infty, |s|]$. Within this interval one can thus replace, for large negative $s$, the Airy function by its asymptotic form for large negative argument $Ai(z) \approx \frac{1}{\sqrt{\pi}}|z|^{-1/4}\sin(\frac{2}{3}|z|^{3/2} + \pi/4)$ for $z \to -\infty$. By substituting the Airy function by its asymptotic behavior in the integral over $u$ in Eq. (126), one finds straightforwardly that $F_{1,b}(s)$ behaves asymptotically as

$$ F_{1,b}(s) \approx \frac{1}{\pi} \sqrt{|s|} \quad \text{for } s \to -\infty, \quad (132) $$

which matches, as it should, with the Wigner semicircle expanded close to $x_{\text{edge}} = \sqrt{2N}/\alpha$ (28).

## 2. Kernel at the edge

We now consider the finite-temperature kernel with both $x$ and $x'$ close to the edge $x_{\text{edge}} = \sqrt{2N}/\alpha$. We set $x = x_{\text{edge}} + w_N s$ and $x' = x_{\text{edge}} + w_N s'$. We insert these coordinates into Eq. (94) and follow the same analysis as in the case of the density (see above). This finally gives in the scaling limit

$$ K_{\mu}(x, x') \approx \frac{1}{w_N} K_b^{\text{edge}} \left( \frac{x - x_{\text{edge}}}{w_N}, \frac{x' - x_{\text{edge}}}{w_N} \right), \quad (133) $$

where the scaled finite-temperature edge kernel is given by

$$ K_b^{\text{edge}}(s,s') = \int_{-\infty}^{\infty} \frac{\text{Ai}(s+u)\text{Ai}(s'+u)}{e^{-bu}+1} du. \quad (134) $$

<figure><img src="image_5.png" /><figcaption>FIG. 4. Plot of the two-point correlation function at the edge $g_b^{\text{edge}}(s, s' = -4)$ as a function of $s$ and for scaled inverse temperatures $b = 0.5$ and 5. Inset: plot of the two-point correlation function in the bulk $g_y^{\text{bulk}}(v)$ versus $v$ for scaled inverse temperatures $y = 0.3, 3$, and 10.</figcaption></figure>

For $s = s'$ it reduces to the expression (126) for the density. Note that in the limit of zero temperature, when $b \to \infty$, the nonzero contribution to the integral over $u$ on the right-hand side of Eq. (134) comes from $u \in [0, +\infty)$ and one gets, using Eq. (37),

$$ \lim_{b \to \infty} K_b^{\text{edge}}(s, s') = \int_0^\infty \text{Ai}(s+u)\text{Ai}(s'+u) du = K_{\text{Airy}}(s, s'). \quad (135) $$

The kernel in Eq. (134) is thus the finite-temperature generalization of the Airy kernel (37). Finally, in Fig. 4 we show a plot of the two-point correlation function at the edge $g_b^{\text{edge}}(s, s') = K_b^{\text{edge}}(s, s)K_b^{\text{edge}}(s', s') - [K_b^{\text{edge}}(s, s')]^2$ for different scaled inverse temperatures $b$.

## C. Extremal statistics near the edge at finite temperature

### 1. Statistics of the rightmost fermion: Exact distribution and its tails

We are now in position to study the fluctuations of the position $x_{\text{max}}(T)$ of the rightmost fermion at finite temperature $T$. In principle, it can be derived from the joint PDF of the fermion positions in Eq. (63) as

$$ \text{Pr}[x_{\text{max}}(T) \le w] = \int_{-\infty}^{w} dx_1 \dots \int_{-\infty}^{w} dx_N P_{\text{joint}}(x_1, \dots, x_N). \quad (136) $$

At $T = 0$, this distribution has a limiting scaling form given by the Tracy-Widom distribution [see Eqs. (47) and (48)], as discussed in Sec. III C. As we have shown in Sec. IV C, in the limit of large $N$, the positions of the fermions at finite $T$ form a determinantal point process, with the kernel parametrized by temperature as given in Eq. (133). As a result (following the discussion in Sec. III C), the multiple integral in Eq. (136) can
be written as a Fredholm determinant [15]. Indeed, at finite temperature $T \sim O(N^{1/3})$ (with $b = \hbar\omega N^{1/3}/T$ fixed) and in the $N \to \infty$ limit, the scaled cumulative distribution function (CDF), denoted by $Q_b(s)$, of $x_{\max}(T)$ can be expressed as the following Fredholm determinant [54]:

$$ \mathrm{Pr} \left[ x_{\max}(T) \le \frac{\sqrt{2N}}{\alpha} + \frac{N^{-\frac{1}{6}}}{\alpha\sqrt{2}} s \right] \underset{N \to \infty}{\to} Q_b(s) \\ := \mathrm{Det}(I - P_s K_b^{\mathrm{edge}} P_s), \qquad (137) $$

where $P_s$ is the projector on the interval $[s, +\infty)$, the kernel $K_b^{\mathrm{edge}}$ is given in Eq. (134), and $\alpha = \sqrt{m\omega/\hbar}$.

The first property to note is that in the limit $T \to 0$, i.e., $b \to +\infty$, since $K_b^{\mathrm{edge}} \to K_{\mathrm{Airy}}$ this distribution (137) converges to the TW distribution for GUE given in Eq. (48). Hence, $Q_b(s)$ is a generalization of the TW distribution to finite temperature. The calculation of the Fredholm determinant (FD) in Eq. (137) is quite involved. As we have pointed out in Ref. [11], the same FD occurs in the exact solution of the KPZ equation with droplet initial conditions. This correspondence is recalled in Sec. VIII B and here we will borrow some of the results obtained in that context. This FD can be expressed in terms of the solution of a nonlocal generalization of the Painlevé II equation [32], namely, one has

$$ \partial_s^2 \ln Q_b(s) = - \int_{-\infty}^{+\infty} dv \sigma_b'(v) [q_b(s, v)]^2, \qquad (138) $$

where

$$ \sigma_b(v) = \frac{1}{1 + e^{-bv}} \qquad (139) $$

and $\sigma_b'(v) = \partial_v \sigma_b(v)$. The function $q_b(s, v)$ satisfies a nonlinear integrodifferential equation in the variable $s$:

$$ \partial_s^2 q_b(s, v) = \left\{ s + v + 2 \int_{-\infty}^{+\infty} dw \sigma_b'(w) [q_b(s, w)]^2 \right\} q_b(s, v), \qquad (140) $$

with the boundary condition $q_b(s, v) \underset{s \to +\infty}{\sim} \mathrm{Ai}(s + v)$. Note that in the zero-temperature limit $b \to \infty$, $\sigma_b'(v) \to \delta(v)$, hence one recovers that $q_b(s, 0)$ satisfies the standard Painlevé II equation (49) which is related to the TW distribution for GUE. The analysis of this equation (140) is rather nontrivial. Alternatively, a numerical evaluation of the FD is possible, along the lines of Ref. [71] using the method developed by Borneman [72]. This is left for future studies.

Tails of $Q_b(s)$. While we have a formal expression for the full scaled distribution function $Q_b(s)$ in terms of a FD in Eq. (137), can we determine its tails for large $|s|$ explicitly for fixed $b$? Indeed, this is possible for arbitrary $b$ for the right tail $s \to +\infty$. However, for the left tail, we can only provide results in the scaling limit $s \to -\infty$, $b \to \infty$, but keeping the ratio $s/b^2$ fixed.

We start with the right tail. In the limit of large positive $s$, the FD in Eq. (137) can be approximated by the first term in a trace expansion. Indeed, in this limit, since $K_b^{\mathrm{edge}}(s, s')$ is “small,” one can expand the FD in Eq. (137) in the following way:

$$ \ln Q_b(s) := \ln \mathrm{Det}(I - P_s K_b^{\mathrm{edge}} P_s) \\ = - \sum_{p=1}^{\infty} \frac{1}{p} \mathrm{Tr} [P_s K_b^{\mathrm{edge}} P_s]^p. \qquad (141) $$

Keeping only the leading $p=1$ term gives, for large $s$,

$$ \ln Q_b(s) \approx -\mathrm{Tr}[P_s K_b^{\mathrm{edge}} P_s] = -\int_s^{+\infty} K_b^{\mathrm{edge}}(s', s') ds'. \qquad (142) $$

Taking derivative w.r.t. $s$, and using $\lim_{s \to \infty} Q_b(s) = 1$, yields

$$ Q_b'(s) \approx K_b^{\mathrm{edge}}(s, s) = F_{1,b}(s) = \int_{-\infty}^{+\infty} du \frac{\mathrm{Ai}(s+u)^2}{1+e^{-bu}}, \qquad (143) $$

where $F_{1,b}(s)$ is the density near the edge in Eq. (126). Therefore, the right tail of $Q_b(s)$ coincides to leading order with the edge density. It turns out that, just the leading term already provides a numerically accurate estimation of the right tail of $Q_b(s)$. Indeed, this is also the case at $T=0$, where the edge density provides a numerically accurate approximation of the right tail of the TW distribution. From the analysis of $F_{1,b}(s)$ in Eq. (131), we see that for fixed $b$, as $s \to \infty$,

$$ Q_b(s) \sim \frac{1}{\sqrt{4\pi b}} \exp\left(-bs + \frac{b^3}{12}\right). \qquad (144) $$

However, if one scales $b$ and $s$ such that $\tilde{s} = s/b^2$ is fixed, then the right tail of $Q_b(s)$ undergoes the same crossover as $F_{1,b}(s)$ at $\tilde{s}_c = \frac{1}{4}$ [as in Eq. (131)].

We now turn to the left tail of $Q_b(s)$ as $s \to -\infty$. Here, analyzing the FD for fixed $b$ with $s \to -\infty$ turns out to be difficult. However, one can make progress in the scaling limit when $b \to \infty$, $s \to -\infty$, keeping the ratio $\tilde{s} = s/b^2$ fixed. In fact, in the context of the height distribution of the KPZ equation at late times (and will be discussed later in Sec. VIII B), this scaling limit was already investigated in Ref. [70] by analyzing the solution of the nonlocal Painlevé equation (140). There it was argued that in this scaling limit the CDF behaves as

$$ Q_b(s) \sim e^{-b^6 \Phi_{-}(s/b^2)}, \quad \text{where} \quad \Phi_{-}(z) = \frac{1}{12} |z|^3. \qquad (145) $$

## 2. Statistics of the rightmost fermion: Finite-temperature behavior of the distribution of the position of the rightmost fermion

In the previous subsection we discussed the limiting distribution of the position of the rightmost fermion in the limit where $T \sim O(N^{1/3})$ and $N$ is large. In that analysis, we kept the scaling parameter $b = \hbar\omega N^{1/3}/T$ fixed and investigated the CDF $Q_b(s)$ as a function of $s$ for fixed $b$. We were able to obtain explicit results for the tails of $Q_b(s)$ for $b$ large, i.e., $T \ll O(N^{1/3})$. Thus, in some sense, the system was still in the vicinity of the $T=0$ limit. In this section, we consider the opposite high-temperature limit where $T \gg O(N^{1/3})$, i.e., the $b \to 0$ limit.

To proceed, we use some recent results from the connection between the KPZ equation in droplet geometry and the fermion problem [11] (for details, see Sec. VIII B). The
height distribution in the KPZ problem was recently analyzed exactly in the short time limit [73], which corresponds to high temperature in the fermion problem. This allows one to obtain the high-temperature expansion, in powers of the small parameter b, of the cumulants of the position of the rightmost fermion $x_{\text{max}}(T)$, and to provide approximate interpolation formula which should be useful for comparison with cold-atom experiments.

*Cumulants.* We start by discussing the mean position and the variance, as obtained from the analysis of Sec. VIII B. Let us define the rescaled variable

$$ \xi = \frac{x_{\text{max}}(T) - x_{\text{edge}}}{w_N}, \quad (146) $$

where $x_{\text{edge}} = \sqrt{2N}/\alpha$ is the edge at $T=0$. Note that this definition differs from the one of the variable $\xi$ defined in Ref. [73]. The first few terms in the series expansion of the mean position read as

$$ \langle \xi \rangle = -\frac{1}{2b} \ln(4\pi b^3) + \frac{\gamma_E}{b} - \sqrt{\frac{\pi}{2}} \frac{b^{1/2}}{2} + \left( \frac{32\pi}{9\sqrt{3}} - 2 - \frac{3\pi}{2} \right) \frac{b^2}{4} + O(b^{7/2}) \quad (147) $$

and the variance behaves as

$$ \begin{aligned} \langle \xi^2 \rangle^c &= \langle \xi^2 \rangle - \langle \xi \rangle^2 \\ &= \frac{\pi^2}{6b^2} + \sqrt{2\pi} \frac{1}{2b^{1/2}} \\ &\quad + \left( 4 + 5\pi - \frac{32\pi}{3\sqrt{3}} \right) \frac{b}{4} + O(b^{5/2}). \end{aligned} \quad (148) $$

These formulas give the behavior for small b, and we know that they should crossover at large b to the zero-temperature limit given by the cumulants of the GUE Tracy-Widom distribution. In particular, the mean $\langle \xi \rangle$ converges to $m_{\text{TW}} = -1.771086807411...$, and the variance $\langle \xi^2 \rangle^c$ converges to $\sigma_{\text{TW}}^2 = 0.8131947928329...$. In fact, we know a bit more: as shown in the Appendix B the large-$b$ expansion takes the form

$$ \langle \xi \rangle = m_{\text{TW}} + O(b^{-4}), \quad (149) $$

$$ \langle \xi^2 \rangle^c = \sigma_{\text{TW}}^2 + \frac{\pi^2}{3b^2} + O(b^{-4}). \quad (150) $$

Note that the form of the correction to the variance as $a_2/b^2$ is common to the KPZ universality class [74], although the prefactor $a_2$ may be nonuniversal. For the continuum KPZ equation its value is fixed and is computed in Appendix B.

Since the variance is the easiest cumulant to measure in a cold-atom experiment, we give here an approximation to the crossover from low to high values of b by constructing a Padé approximant $\langle \xi^2 \rangle^c|_{\text{Padé}}$ which (i) reproduces all known terms in the small-$b$ expansion given in Eq. (148) and (ii) converges at large b to $\sigma_{\text{TW}}^2$ as $\langle \xi^2 \rangle^c|_{\text{Padé}} \simeq \sigma_{\text{TW}}^2 + \frac{a_2}{b^2}$. It reads as, with

<figure><img src="image_6.png" /><figcaption>FIG. 5. Plot of the Padé approximate $\langle \xi^2 \rangle^c|_{\text{Padé}}$ given in Eq. (151) with $a_2 = \pi^2/3$. The horizontal line corresponds to the variance of the TW distribution $\sigma_{\text{TW}}^2$, which is the exact result for $b \to \infty$, i.e., in the $T \to 0$ limit.</figcaption></figure>

the value $a_2 = \pi^2/3$,

$$ \langle \xi^2 \rangle^c|_{\text{Padé}} = \frac{1.64493 + 1.13494 b^{3/2} + a_2 b^4 + 0.813195 b^6}{b^2(1 - 0.0719632 b^{3/2} + b^4)} \quad (151) $$

and is plotted in Fig. 5. Although there is some degree of arbitrariness, this curve should be useful to calibrate the experiments, given that the range of values of b presently available is $b \sim 0.5-2$. Another way to present the result for the variance of the fluctuations of the rightmost fermion is to divide it by the (half) size of the Fermi cloud (at $T=0$) $x_{\text{edge}}$, and write

$$ \left[ \frac{x_{\text{max}}(T) - x_{\text{edge}}}{x_{\text{edge}}} \right]^2 \right]^{c} = \frac{T^2}{T_F^2} \mathcal{V} \left( b = \frac{T_F}{N^{2/3} T} \right), \quad (152) $$

where $T_F = N\hbar\omega$, and the scaling function $\mathcal{V}(z)$ is such that $\mathcal{V}(0) = \frac{\pi^2}{24}$ [see Eq. (148)] and $\mathcal{V}(z) \simeq z \to +\infty \sigma_{\text{TW}}^2 z^2/4$ to yield back the zero-temperature limit (i.e., as $b \to \infty$). A Padé approximation of the function $\mathcal{V}(z)$ is easily obtained from (151).

Finally, we also give the third cumulant of the scaled variable $\xi$ in Eq. (281). For small b it reads as

$$ \langle \xi^3 \rangle^c = \frac{2\zeta(3)}{b^3} + \left( \frac{32}{3\sqrt{3}} - 6 \right) \frac{\pi}{4} + O(b^{3/2}) \quad (153) $$

which allows to calculate the skewness of the position of the rightmost fermion (note that the skewness is independent of any rescaling)

$$ \begin{aligned} \text{Sk} := \frac{\langle \xi^3 \rangle^c}{[\langle \xi^2 \rangle^c]^{3/2}} &= 1.13955 - 1.30237 b^{3/2} \\ &\quad + 1.20563 b^3 + O(b^{7/2}). \end{aligned} \quad (154) $$

It decreases from the skewness of the Gumbel distribution (see below) $\text{Sk}_{\text{Gumbel}} = 1.13955$ at high temperature (small b) to the skewness of the TW distribution for GUE $\text{Sk}_{\text{TW}} = 0.224084203610...$ at zero temperature (large b). Note, however, that since $\langle \xi^3 \rangle^c = \kappa_3 + O(b^{-4})$ (see Appendix B),
we find that the approach to the zero-temperature limit is from below

$$ Sk = Sk_{TW} \left[ 1 - \frac{\pi^2}{2\sigma_{TW}^2 b^2} + O(b^{-4}) \right]. \quad (155) $$

Hence, the skewness is (weakly) nonmonotonic as a function of the temperature.

**High-temperature limit and the Gumbel distribution.** In the high-temperature limit of the edge regime, i.e., $T \gg N^{1/3} \hbar\omega$ (equivalently $b \ll 1$), the full PDF of $\xi$ becomes a Gumbel distribution, up to a temperature-dependent constant shift [65,73]. More precisely, one finds

$$ \frac{x_{\max}(T) - x_{\text{edge}}}{x_{\text{edge}}} \approx \frac{1}{4y} \ln \left( \frac{N^2}{4\pi y^3} \right) + \frac{1}{2y} \gamma, \quad y = \frac{\hbar\omega N}{T} = \frac{T_F}{T}, \quad (156) $$

where $\gamma$ is a Gumbel random variable with PDF $p(\gamma) = e^{-\gamma}e^{-\gamma}$. Note that we have used the variable $y$ which was introduced to study the high-temperature bulk regime [which is defined by $y = O(1)$]. Although this formula is derived here from studying the high-temperature limit of the low-temperature edge regime, we expect that it should match in some way with the low-temperature limit of the high-temperature bulk regime, although as we will see this matching is not trivial. This result for high temperature in Eq. (156) together with the $T=0$ result [see Eqs. (47) and (48)] show that the fluctuations of the position of the rightmost fermion, in the edge regime, interpolates between a TW-random variable $\chi_2$ when $T \ll N^{1/3} \hbar\omega$ and Gumbel random variable $\gamma$ when $T \gg N^{1/3} \hbar\omega$ [65,73], i.e. (in terms of the original parameters $T$ and $N$),

$$ \frac{x_{\max}(T) - x_{\text{edge}}}{x_{\text{edge}}} \sim \begin{cases} \frac{1}{2} N^{-2/3} \chi_2, & T \ll N^{1/3} \hbar \omega, \\ \frac{T}{4N \hbar \omega} \ln \left[ \frac{T^3}{4\pi N(\hbar \omega)^3} \right] + \frac{T}{2N \hbar \omega} \gamma, & T \gg N^{1/3} \hbar \omega. \end{cases} \quad (157) $$

It is instructive to compare this Gumbel distribution to the one that we can infer, from extreme value statistics arguments (see e.g., Ref. [75]), in the very high-temperature regime $y \ll 1$. In that regime, we know that the positions of the fermions are independent and identical random variables drawn from the Boltzmann distribution given in Eq. (98), and therefore $x_{\max}(T)$ is also distributed according to a Gumbel law. This implies

$$ \frac{x_{\max}(T) - x_{\text{edge}}}{x_{\text{edge}}} \approx \sqrt{\frac{\ln N}{y}} - 1 - \frac{\ln(4\pi \ln N)}{4\sqrt{y \ln N}} + \frac{1}{2\sqrt{y \ln N}} \tilde{y}, \quad (158) $$

where $\tilde{y}$ is also a Gumbel random variable [we used the different notation $\tilde{y}$ in Eq. (158) to emphasize that this random variable is different from the Gumbel variable in Eq. (156), while both of them have the same statistics]. The behavior (158) is expected to hold for $N \gg 1$, and fixed $y \ll 1$. Note that $x_{\text{edge}}$ in this formula is the $T=0$ edge, while at these

<figure><img src="image_7.png" /><figcaption>FIG. 6. Sketch of the behavior of the distribution of the rightmost fermion as a function of $T/\hbar\omega$. When $T/\hbar\omega \ll N^{1/3}$, the distribution is given by the TW distribution for GUE [see Eq. (157)] while for $T/\hbar\omega \gg N^{1/3}$, and $T/\hbar\omega \ll N$, it crosses over to a Gumbel distribution (157), denoted as Gumbel (I) in the figure. The full crossover, when $T/\hbar\omega \sim O(N^{1/3})$ is described the finite-temperature generalization of the TW distribution (137), which also appears in the KPZ equation in a droplet geometry [see also Eq. (151) for the second cumulant]. Finally, for $T/\hbar\omega \gg N$, the distribution is described by yet another Gumbel distribution [denoted as Gumbel (II) in the figure], as discussed in the text [see Eq. (158)].</figcaption></figure>

temperatures the density does not display a true edge, and extends well beyond $x_{\text{edge}}$. We note that, although we obtain Gumbel laws in the two regimes (156) and (158), the detailed dependence in $y$ and $N$ is quite different. From its derivation, the regime (156) should hold for large $y$. If one compares the deterministic terms in Eqs. (156) and (158), keeping in mind that the fluctuations increase with temperature, one obtains the stronger condition for (156) to hold:

$$ y \gg \ln N, \quad \text{i.e.,} \quad N^{1/3} \ll \frac{T}{\hbar\omega} \ll \frac{N}{\ln N}. \quad (159) $$

The interpolation between these two regimes remains an open problem. It requires to study the Fredholm determinant associated with the full finite-temperature kernel (94) in the region $y = O(1)$. Note that, although the fluctuations of $x_{\max}(T)$ are universal (this is shown in Sec. VII D) in the low-temperature scaling regime $b = O(1)$ ([i.e., in Eq. (157)], the extreme value statistics arguments leading to (158) also show that the (logarithmic) dependence in $N$ depends on the power-law index $p$ of the confining potential $V(x) \sim |x|^p$. These different behaviors of the rightmost fermion, as described in Eqs. (157) and (158) above, are sketched in Fig. 6.

Note that aside from the regime of typical fluctuations of $x_{\max}(T)$ one can also study the large deviations. In the limit of high temperature of the low-temperature edge regime, this was done in detail in Ref. [73]. In addition, one can also study the counting statistics of the number of fermions in a given interval $J$, or the PDF of the spacing between nearest-neighbor fermions near the edge at finite temperature. As a consequence of the determinantal structure of all correlations in the large-$N$ limit, these observables are given in terms of Fredholm determinants by exactly the same formula as (39) and (42) replacing the kernel $K_J$ there by the finite-$T$ edge kernel (134). The detailed analysis of such FD is left for future investigations.

Here, we add a few interesting remarks. First one can extend the property shown in Eq. (54) to finite temperature, assuming that again the equivalence canonical-grand canonical goes through at large $N$ (the property is exact in the grand-canonical ensemble). It can be written everywhere (bulk and edge) but let us display it here near the edge. Defining the rescaled positions of the fermions near the edge $a_i = \frac{x_i - x_{\text{edge}}}{w_N}$ and taking
the $N \to \infty$ limit, keeping $a_i$'s fixed

$$ 
\begin{aligned}
\left\langle \prod_{i=1}^{\infty} f(a_i) \right\rangle_T &= \mathrm{Det}(I - L_{f,b}), \quad L_{f,b}(r,r') \\
&= [1 - f(r)] K_b^{\mathrm{edge}}(r,r'),
\end{aligned}
\quad (160) $$

where $K_b^{\mathrm{edge}}(r,r')$ is given in Eq. (134). In Eq. (160), the average is performed at finite temperature $T$ and the function $f(x)$ is arbitrary, provided the right-hand side exists. Next, since one can rewrite the CDF of the position of the rightmost fermion in Eq. (137) more explicitly in terms of the Airy kernel as

$$ Q_b(s) := \mathrm{Det}(I - P_s K_b^{\mathrm{edge}} P_s) = \mathrm{Det}(I - \bar{K}_{b,s}), \quad (161) $$

$$ \bar{K}_{b,s}(u, u') = K_{Ai}(u, u') \sigma_{s,b}(u), \quad \sigma_{s,b}(u) = 1/(1+e^{-b(u-s)}), \quad (162) $$

we immediately obtain [with $P_J(x)$ denoting the indicator function of the interval $J$]

$$ 
\begin{aligned}
\left\langle \prod_{i=0}^{\infty} P_{(-\infty, s]}(a_i) \right\rangle_T &= Q_b(s) = \left\langle \prod_{i=0}^{\infty} [1 - \sigma_{s,b}(a_i)] \right\rangle_{T=0} \\
&= \left\langle \prod_{i=0}^{\infty} \frac{1}{1 + e^{b(a_i - s)}} \right\rangle_{T=0}.
\end{aligned}
\quad (163) $$

A similar product appeared in a recent paper by Borodin and Gorin [76].

## VI. THE D-DIMENSIONAL ISOTROPIC HARMONIC OSCILLATOR AT ZERO TEMPERATURE

We now consider the model for noninteracting fermions in a harmonic potential in arbitrary dimension $d$, as defined in Sec. I B. Here, we focus on $T = 0$, the finite-$T$ case is studied in the next section. To study this more general problem it is convenient to use a method based on the one-body Euclidean propagator that we describe in Sec. V A. This will lead to results both in the bulk (Sec. V B), as well as at the edge (Sec. V C) of the $d$-dimensional Fermi gas. In particular, this provides a useful alternative (also in $d = 1$) to the method relying on the Plancherel-Rotach asymptotic formulas for Hermite polynomials.

### A. Representation of the $T = 0$ kernel using the one-body Euclidean propagator in arbitrary $d$

#### 1. General framework: Kernel and propagator

We start by considering the zero-temperature kernel for $N$ noninteracting spinless fermions with an arbitrary one-body Hamiltonian $\hat{H}$, as defined in Eq. (1). As discussed in Sec. II, the kernel corresponding to a system with Fermi energy $\mu$ is then defined 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 (164) $$

where $\theta(z)$ is the Heaviside step function and the energy eigenvalues $\epsilon_k$ are labeled by $d$ quantum numbers denoted

by $\mathbf{k} \in \mathbb{Z}^d$. We now compute

$$ 
\begin{aligned}
G(\mathbf{x}, \mathbf{y}; t) &= \frac{t}{\hbar} \int_0^\infty d\mu \exp\left(-\frac{t\mu}{\hbar}\right) K_\mu(\mathbf{x}, \mathbf{y}) \\
&= \sum_{\mathbf{k}} \psi_{\mathbf{k}}^*(\mathbf{x}) \psi_{\mathbf{k}}(\mathbf{y}) \exp\left(-\frac{\epsilon_{\mathbf{k}}t}{\hbar}\right)
\end{aligned}
\quad (165) $$

and immediately see that $G(\mathbf{x}, \mathbf{y}; t)$ is in fact the one-body Euclidean propagator associated to the one-body Hamiltonian $\hat{H}$. By definition, it obeys the imaginary-time Schrödinger equation

$$ -\hbar \frac{\partial G(\mathbf{x}, \mathbf{y}; t)}{\partial t} = \hat{H} G(\mathbf{x}, \mathbf{y}; t), \quad (166) $$

where $\hat{H} = \hat{H}(\mathbf{y}, \frac{\hbar}{i}\nabla_{\mathbf{y}})$ is the quantum Hamiltonian defined in Eq. (1) acting on the variable $\mathbf{y}$ (in our convention). From the completeness of the basis of the eigenfunctions, it satisfies the initial condition

$$ G(\mathbf{x}, \mathbf{y}; 0) = \delta^d(\mathbf{x} - \mathbf{y}). \quad (167) $$

If the propagator $G$ is known as a function of $t$, then the kernel can be obtained via the Bromwich inversion formula for Laplace transforms as

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

where $\Gamma$ indicates the Bromwich integration contour in the complex plane.

For the isotropic $d$-dimensional harmonic oscillator, which we focus on here, the Euclidean propagator is known exactly at all times, and is given by [77]

$$ G(\mathbf{x}, \mathbf{y}; t) = \left[ \frac{\alpha^2}{2\pi \sinh(\omega t)} \right]^{d/2} \exp \left\{ -\frac{\alpha^2}{2\sinh(\omega t)} \times \left[ (\mathbf{x}^2 + \mathbf{y}^2)\cosh(\omega t) - 2\mathbf{x} \cdot \mathbf{y} \right] \right\}, \quad (169) $$

where $\alpha = \sqrt{m\omega/\hbar}$. The Bromwich integral in Eq. (168) is in general difficult to compute explicitly. However, for a system with a large number of fermions, the Fermi energy $\mu$ will be large and so we can assume that the integral is dominated by small values of $t$. The validity of the approach will be made more precise below.

#### 2. Short-time expansion of the propagator and of the kernel

As we will make clear later, to obtain all of our results, we will need the short-time expansion of the propagator up to $O(t^3)$. For this purpose it is useful to write the propagator $G(\mathbf{x}, \mathbf{y}; t)$ given in Eq. (169) in the following form:

$$ G(\mathbf{x}, \mathbf{y}; t) = \left[ \frac{\alpha^2}{2\pi \sinh(\omega t)} \right]^{d/2} \exp \left( -\frac{\alpha^2}{2\sinh(\omega t)} \times \left\{ (\mathbf{x} - \mathbf{y})^2 + (\mathbf{x}^2 + \mathbf{y}^2)[\cosh(\omega t) - 1] \right\} \right), \quad (170) $$
which gives at coinciding points with $r = |\mathbf{x}|$

$$ G(\mathbf{x}, \mathbf{x}; t) = \left[ \frac{\alpha^2}{2\pi \sinh(\omega t)} \right]^{\frac{d}{2}} \exp \left[ -\tanh \left( \frac{\omega t}{2} \right) \alpha^2 r^2 \right]. \quad (171) $$

Expanding the formula (170) to $O(t^3)$ gives

$$ G(\mathbf{x}, \mathbf{y}; t) \approx \left(\frac{m}{2\pi\hbar t}\right)^{d/2} \exp \left\{ -\frac{m}{2\pi\hbar t} (\mathbf{x}-\mathbf{y})^2 - \frac{m\omega^2 t}{12\hbar} [3(\mathbf{x}^2 + \mathbf{y}^2) - (\mathbf{x}-\mathbf{y})^2] - \frac{d\omega^2 t^2}{12} + \frac{m\omega^4 t^3}{720\hbar} [15(\mathbf{x}^2 + \mathbf{y}^2) - 7(\mathbf{x}-\mathbf{y})^2] \right\}. \quad (172) $$

In particular, at coinciding points $\mathbf{y} = \mathbf{x}$ with $r = |\mathbf{x}|$ and using $\alpha = \sqrt{m\omega/\hbar}$, one obtains to $O(t^3)$

$$ G(\mathbf{x}, \mathbf{x}; t) \approx \left[ \frac{\alpha^2}{2\pi\omega t} \right]^{\frac{d}{2}} \exp \left( -\frac{\alpha^2 r^2 \omega t}{2} - \frac{d\omega^2 t^2}{12} + \frac{\alpha^2 r^2 \omega^3 t^3}{24} \right). \quad (173) $$

Substituting the expression (173) into Eq. (168) then yields the short-time expansion to order $O(t^3)$ of the kernel:

$$ K_{\mu}(\mathbf{x}, \mathbf{y}) \approx \left[ \frac{\alpha^2}{\omega} \right]^{\frac{d}{2}} \int_{\Gamma} \frac{dt}{(2\pi t)^{1+\frac{d}{2}i}} \exp \left( -\frac{m}{2\pi\hbar t}(\mathbf{x}-\mathbf{y})^2 + \frac{t}{\hbar} \left\{ \mu - \frac{m\omega^2}{12}[3(\mathbf{x}^2+\mathbf{y}^2) - (\mathbf{x}-\mathbf{y})^2] \right\} \right. \\ \left. - \frac{d\omega^2 t^2}{12} + \frac{m\omega^4 t^3}{720\hbar}[15(\mathbf{x}^2+\mathbf{y}^2) - 7(\mathbf{x}-\mathbf{y})^2] \right). \quad (174) $$

At coinciding points, the kernel expanded to $O(t^3)$ becomes

$$ K_{\mu}(\mathbf{x}, \mathbf{y}) \approx \left[ \frac{\alpha^2}{\omega} \right]^{\frac{d}{2}} \int_{\Gamma} \frac{dt}{(2\pi t)^{1+\frac{d}{2}i}} \exp \left( \frac{t}{\hbar} \left[ \mu - \frac{m\omega^2 r^2}{2} \right] - \frac{d\omega^2 t^2}{12} + \frac{\alpha^2 r^2 \omega^3 t^3}{24} \right). \quad (175) $$

The latter two formulas will be used extensively below.

## B. Results in the bulk

### 1. Bulk density profile

We start by analyzing the kernel at coinciding points, i.e., $\mathbf{x} = \mathbf{y}$, to obtain the density $\rho_N(\mathbf{x}) = (1/N)K_\mu(\mathbf{x}, \mathbf{x})$ for large $N$ (equivalently for large $\mu$ as we assume now and verify a *posteriori*). Our starting point is Eq. (175) where we see that the effective energy scale entering the Laplace transform is not simply $\mu$ but

$$ \epsilon(r) = \mu - m\omega^2 r^2 / 2 = \mu \left( 1 - \frac{r^2}{r_{\text{edge}}^2} \right), \quad (176) $$

which is the difference between the Fermi energy and the classical potential energy at the radial coordinate $r$. We have defined

$$ r_{\text{edge}} = \sqrt{2\mu/m\omega^2}, \quad (177) $$

the radius at which $\epsilon(r)$ vanishes. We will see below that $r_{\text{edge}}$ is also the edge where the bulk density vanishes. When $\epsilon(r)$ is large, which is the case for large $\mu$ and away from the edge, the integral in Eq. (175) can be approximated by the term of $O(t)$ in the exponential, and is dominated by times of order $t \sim \hbar/\epsilon(r)$. The corrections to this approximation coming from the $O(t^2)$ and $O(t^3)$ are, respectively,

$$ \frac{d}{12} \left[ \frac{\hbar\omega}{\epsilon(r)} \right]^2 + \frac{\alpha^2 r^2}{24} \left[ \frac{\hbar\omega}{\epsilon(r)} \right]^3. \quad (178) $$

Hence, the first one is unimportant when

$$ \epsilon(r) \gg \hbar\omega \iff |r - r_{\text{edge}}| \gg (r_{\text{edge}}\alpha)^{-1/\alpha}. \quad (179) $$

Similarly, the cubic terms can be seen to be unimportant when

$$ \frac{\alpha^2 r^2}{24} \left[ \frac{\hbar\omega}{\epsilon(r)} \right]^3 \ll 1 \iff |r - r_{\text{edge}}| \gg (r_{\text{edge}}\alpha)^{-1/\alpha}. \quad (180) $$

We will see below that $r_{\text{edge}}$ is large for large $N$. In that case, sufficiently away from the edge, Eq. (180) is satisfied and this automatically ensures that Eq. (179) is also satisfied. Hence, keeping only this leading term to describe the bulk density, we can use the standard identity

$$ \int_{\Gamma} \frac{dt}{2\pi i} \frac{1}{t^{d/2+1}} \exp(zt) = \frac{z^{d/2}}{\Gamma(1+\frac{d}{2})} \theta(z) \quad (181) $$

to write

$$ K_{\mu}(\mathbf{x}, \mathbf{y}) \approx \frac{1}{\Gamma(1+\frac{d}{2})} \left[ \frac{m}{2\pi\hbar^2} \right]^{\frac{d}{2}} \left( \mu - \frac{m\omega^2 r^2}{2} \right)^{\frac{d}{2}} \theta \left( \mu - \frac{m\omega^2 r^2}{2} \right) \quad (182) $$

$$ = \frac{1}{2^d \pi^{\frac{d}{2}} \Gamma(1 + \frac{d}{2})} \alpha^{2d} (r_{\text{edge}}^2 - r^2)^{\frac{d}{2}} \theta(r_{\text{edge}} - r), \quad (183) $$

where $\theta$ is the Heaviside step function. Consequently, the bulk density, as a function of $r = |\mathbf{x}|$ is given by

$$ \rho_N(r) = \frac{1}{2^d N \pi^{\frac{d}{2}} \Gamma(1 + \frac{d}{2})} \alpha^{2d} (r_{\text{edge}}^2 - r^2)^{\frac{d}{2}} \theta(r_{\text{edge}} - r). \quad (184) $$

Hence, as anticipated above, the bulk density vanishes at the radius $r = r_{\text{edge}}$ given in Eq. (177). It should also be noticed
that the expression given here for the density is the one obtained by the local density approximation (LDA) [4].

The value of the chemical potential $\mu$ corresponding to a fixed value of $N$, the number of fermions, can be evaluated from the normalization condition $\int d^d \mathbf{x} \rho_N(\mathbf{x}) = 1$ which yields

$$ \frac{S_d r_{\text{edge}}^d}{\Gamma\left(1 + \frac{d}{2}\right)} \mu^d \left(\frac{m}{2\pi\hbar^2}\right)^{\frac{d}{2}} \int_0^1 du \, u^{d-1} (1-u^2)^{\frac{d}{2}} = N, \quad (185) $$

where $S_d = 2\pi^{d/2}/\Gamma(d/2)$ is the surface unit area of a unit sphere in $d$ dimensions. Using

$$ \int_0^1 du \, u^{d-1} (1-u^2)^{\frac{d}{2}} = \frac{1}{2} B\left(\frac{d}{2}, \frac{d}{2} + 1\right), \quad (186) $$

where $B(p,q)$ denotes the beta function, as well as Eqs. (177) and (185), we obtain the Fermi level and the value of the edge radius as a function of $N$:

$$ \mu = \hbar\omega[N\Gamma(d+1)]^{\frac{1}{d}}, \quad r_{\text{edge}} = \frac{2^{\frac{1}{2}}[\Gamma(d+1)]^{\frac{1}{2d}}}{\alpha} N^{\frac{1}{2d}}, \quad (187) $$

as anticipated in Eq. (6). Finally, this leads to the final result for the normalized density as a function of $N$:

$$ \rho_N(\mathbf{x}) = \frac{1}{2^d N^{\frac{1}{2}} \pi^{\frac{d}{2}} \Gamma\left(1 + \frac{d}{2}\right)} \left(\frac{m\omega}{\hbar}\right)^{\frac{d}{2}} \times \left[ 2\Gamma(d+1)^{\frac{1}{d}} - \frac{m\omega}{\hbar N^{\frac{1}{d}}} r^2 \right]^{\frac{d}{2}}. \quad (188) $$

Using the explicit dependence of $\mu$ and thus $r_{\text{edge}}$ on $N$ from (187) in Eq. (180), we see that this formula for the bulk density is valid in the regime $r_{\text{edge}} - r \gg N^{-1/d}/\alpha$. Below, we will see how to obtain both the density and the kernel in the edge regime where $r_{\text{edge}} - r \sim N^{-1/d}/\alpha$.

As an example, we consider the simple harmonic oscillator in two dimensions. We have numerically computed the density $\rho_N(\mathbf{x})$ for $N = 28$ (corresponding to a full Fermi level) free fermions in a two-dimensional harmonic trap for $\alpha = \sqrt{m\omega/\hbar} = 1$, by directly summing the modulus squared of the first 28 wave functions. Shown in Fig. 7 is the comparison of the numerical summation with the bulk asymptotic formula (188).

## 2. Kernel in the bulk

In order to compute the bulk kernel for the isotropic simple harmonic oscillator, we extend our analysis for the bulk density, and thus we keep only terms up to order $O(t)$ in Eq. (174) which gives

$$ K_{\mu}(\mathbf{x}, \mathbf{y}) \approx \left(\frac{m}{2\pi\hbar}\right)^{\frac{d}{2}} \int_{\Gamma} \frac{dt}{2\pi i} \frac{1}{t^{\frac{d}{2}+1}} \exp\left[\frac{\mu t}{\hbar} - \frac{m(\mathbf{x} - \mathbf{y})^2}{2\hbar t} - \frac{m\omega^2(\mathbf{x}^2 + \mathbf{y}^2 + \mathbf{x} \cdot \mathbf{y})t}{6\hbar}\right]. \quad (189) $$

To proceed, one can use the following integral representation of the standard Bessel function of the first kind of index $\nu$,

<figure><img src="image_8.png" /><figcaption>FIG. 7. Plot of the bulk density (solid line) $\rho_{28}(r)$ for 28 fermions in an isotropic harmonic potential in two dimensions compared with the asymptotic formula for the bulk density (188) (the region below the asymptotic result is shown via the shaded region).</figcaption></figure>

denoted by $J_v(x)$ [78]:

$$ \int_{\Gamma} \frac{dt}{2\pi i} \frac{1}{t^{d/2+1}} \exp\left(zt - \frac{a}{t}\right) = \left(\frac{z}{a}\right)^{d/4} J_{\frac{d}{2}}(2\sqrt{az}). \quad (190) $$

We use this identity to evaluate the representation of the kernel in Eq. (189) to obtain

$$ K_{\mu}(\mathbf{x}, \mathbf{y}) \approx \left(\frac{\alpha^2}{2\pi}\right)^{\frac{d}{2}} \left( \frac{3r_{\text{edge}}^2 - \mathbf{x}^2 - \mathbf{y}^2 - \mathbf{x} \cdot \mathbf{y}}{3(\mathbf{x} - \mathbf{y})^2} \right)^{\frac{d}{4}} \times J_{\frac{d}{2}}\left(\frac{\alpha^2}{\sqrt{3}}|\mathbf{x} - \mathbf{y}|\sqrt{3r_{\text{edge}}^2 - \mathbf{x}^2 - \mathbf{y}^2 - \mathbf{x} \cdot \mathbf{y}}\right), \quad (191) $$

where we recall that $\mu = \frac{1}{2}m\omega^2 r_{\text{edge}}^2$. Now, if we consider two points $\mathbf{x}'$ and $\mathbf{y}'$ both close to a point $\mathbf{x}$ in the bulk we find that the kernel has the scaling form

$$ K_{\mu}(\mathbf{x} + \mathbf{x}', \mathbf{y} + \mathbf{y}') \approx \ell(\mathbf{x})^{-d} K_{d}^{\text{bulk}}[|\mathbf{x}' - \mathbf{y}'|/\ell(\mathbf{x})], \quad (192) $$

where

$$ \ell(\mathbf{x}) = [N\rho_N(\mathbf{x})\gamma_d]^{-1/d}, \quad \gamma_d = S_d/d = \pi^{d/2}[\Gamma(d/2 + 1)] \quad (193) $$

is the local typical separation between fermions in the bulk. The explicit formula for the scaling function in Eq. (192) is given by

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

which has a well-defined limit at the origin with $K_d^{\text{bulk}}(0) = 1/\gamma_d$. In $d=1$, using $J_{1/2}(z) = \sqrt{2/(\pi z)} \sin z$, we recover the standard sine kernel $K_{\text{bulk}}(x) = \frac{\sin(2x)}{\pi x}$ of RMT. The domain of validity of the bulk analysis carried out above can be determined by examining the corrections coming from the quadratic and cubic corrections in $t$, as was done above for the study of the density. The result for the bulk kernel in Eq. (189) is found to be valid when both $\mathbf{x}$ and $\mathbf{y}$ are in the bulk, away from the edge $r_{\text{edge}}$.

The result for the kernel (194) for relative distances of order the typical particle separation can also be obtained from the
LDA [4,79]. Indeed, we remark here that the Fourier transform of this kernel is Fermi-step function, this can be seen using the formula

$$ \int_{|\mathbf{k}|<k_F} \frac{d^d k}{(2\pi)^d} e^{i\mathbf{k}\cdot\mathbf{x}} = \left(\frac{k_F}{2\pi|\mathbf{x}|}\right)^{d/2} J_{d/2}(k_F|\mathbf{x}|). \quad (195) $$

In the expression in Eq. (195) the local Fermi momentum must be chosen as $k_F = k_F(\mathbf{x}) = 2/\ell(\mathbf{x}) = 2[N\rho_N(\mathbf{x})\gamma_d]^{1/d}$. This value is exactly consistent with the one obtained by the counting of states for a uniform system of density $\rho_N(\mathbf{x})$, by setting $N \rho_N(\mathbf{x}) = \int_{|\mathbf{k}|<k_F} \frac{d^d k}{(2\pi)^d} = k_F^d/(4\pi)^{d/2}\Gamma(1+d/2)$. Thus, to describe correlations on scale $\ell(\mathbf{x})$ in the bulk, one can approximate locally the system by free fermions without any external potential, but at a fixed density $\rho_N(\mathbf{x})$, assumed to be slowly varying on that scale. Our method thus gives a rigorous derivation of the results, within the bulk, obtained from the heuristic LDA.

### C. Results at the edge

#### 1. Density near the edge

We have seen how the density and kernel in the bulk region can be obtained via a leading-order short-time expansion of the propagator for the simple harmonic oscillator in Sec. VI B. However, this expansion becomes insufficient when the inequalities in Eqs. (179) and (180) are violated: this occurs near the edge where the density vanishes. As demonstrated below, the description of the edge regime requires one to expand the propagator to higher orders in $t$ as in Eq. (173).

We start by investigating the density near $r = r_{\text{edge}} = \sqrt{2\mu/m\omega^2}$, for finite but large $N$, by setting

$$ r = |\mathbf{x}| = r_{\text{edge}} + z w_N, \quad w_N = b_d N^{-\phi}, \quad (196) $$

where the distance from the edge is parametrized by the scaled dimensionless distance $z$, while the edge exponent $\phi$ is to be determined a posteriori. For convenience, we introduced the factor $b_d = [\Gamma(1+d)]^{-1/d}/(\alpha\sqrt{2})$ which has the dimension of a length.

To calculate the kernel in this edge region, we need to use the expansion of the kernel $K_\mu(\mathbf{x}, \mathbf{x})$ up to $O(t^3)$, as given in Eq. (175). Substituting (196) into (175), we see that the two terms of $O(t)$ cancel each other inside the argument of the exponential, leaving only three terms

$$ K_{\mu}(\mathbf{x}, \mathbf{x}) \approx \left(\frac{\alpha^2}{2\pi\omega}\right)^{d/2} \int_{\Gamma} \frac{dt}{2\pi i} \frac{1}{t^{d/2+1}} \\ \times \exp\left[-\sqrt{\frac{2\mu}{m}}\alpha^2 z b_d N^{-\phi}t - \frac{d}{12}\omega^2 t^2 + \frac{\mu\omega^2}{12\hbar}t^3\right]. \quad (197) $$

We now determine the exponent $\phi$ by comparing the magnitude of the three terms inside the exponential and using $\mu \sim N^{1/d}$. We obtain, in the order that they appear,

$$ T_1 \sim z\sqrt{\mu}N^{-\phi}t \sim zN^{\frac{1}{2d}-\phi}t, \quad T_2 \sim t^2, \quad T_3 \sim \mu t^3 \sim N^{\frac{1}{d}}t^3. \quad (198) $$

Since the first term must be of order $O(1)$, this implies that $t \sim N^{\phi-1/2d}$. We then have only two possibilities a priori. The first

one is to choose $\phi = 1/(2d)$ that keeps $T_2 \sim O(1)$, but then $T_3 \sim N^{1/d}$ and diverges as $N \to \infty$, which is inconsistent. Hence, in order to make the term $T_3 = O(1)$, we must choose

$$ \phi = \frac{1}{6d} \Rightarrow t \sim N^{-1/3d}, \quad (199) $$

which is consistent with the assumption of an expansion in small $t$. It is easy to check in Eq. (171) that with this scaling exponent the terms of order higher than $O(t^3)$ vanish as $N \to \infty$. Note that this result can be understood qualitatively by arguing that there should be of order one particle in the typical fluctuation region, i.e., in a box of linear size $w_N$ around the edge, which leads to

$$ \left(\frac{w_N}{r_{\text{edge}}}\right)^{d-1} \int_{r_{\text{edge}}-w_N}^{r_{\text{edge}}} \rho_N(r) r^{d-1} dr \sim \frac{1}{N} \quad (200) $$

which in turn implies $w_N \sim N^{-1/(6d)}$, using the formula (188) for the density $\rho_N(r) = \rho_N(|\mathbf{x}|)$. Hence, rescaling $t = N^{-1/(3d)}/\{\omega[\Gamma(d+1)]^{1/(3d)}\}\tau$, we obtain our main result for the density

$$ N\rho_N(\mathbf{x}) = K_\mu(\mathbf{x}, \mathbf{x}) = \frac{1}{w_N^d} F_d\left(\frac{r - r_{\text{edge}}}{w_N}\right), \quad (201) $$

where we have defined the width of the edge regime in $d$ dimension

$$ w_N = b_d N^{-1/6d} = \frac{1}{\alpha\sqrt{2}}[\Gamma(1+d)N]^{-1/6d} \quad (202) $$

and the scaling function $F_d(z)$ 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 (203) $$

The expression (203) can be rewritten in the following way. First, we make use of the identity

$$ \frac{1}{\tau^{d/2+1}} = \frac{1}{\Gamma(d/2+1)} \int_0^\infty e^{-\tau x} x^{d/2} dx \quad (204) $$

in Eq. (203) to obtain

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

Rescaling $\tau \to 2^{2/3}\tau$ and using the integral representation of the Airy function

$$ \text{Ai}(z) = \int_{\Gamma} \frac{d\tau}{2\pi i} e^{-\tau z + \tau^3/3} \quad (206) $$

then gives, after a further rescaling $x \to 2^{-2/3}x$, in the integral expression for $F_d(z)$ in Eq. (205) then yields

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

General properties of the density scaling function in $d = 1,2,3$. This scaling function satisfies some general properties in any $d$. For instance, differentiating Eq. (207) w.r.t. $z$ once and using integration by parts, one finds the recursion formula

$$ \frac{dF_d(z)}{dz} = -\frac{1}{4\pi} F_{d-2}(z), \quad (208) $$
which allows to obtain $F_d(z) = \frac{1}{4\pi} \int_z^{+\infty} F_{d-2}(z)$ from the knowledge of $F_{d-2}(z)$. In addition, it is easy to see that $F_d(z)$ satisfies a third-order differential equation

$$ \frac{d^3 F_d(z)}{dz^3} - 4z \frac{d F_d(z)}{dz} + 2d F_d(z) = 0, \quad (209) $$

which must be complemented by appropriate boundary conditions (see below).

Explicit forms of the density scaling function in $d = 1, 2, 3$. In $d = 1$, the above integral can be performed exactly. We start with the identity [80]

$$ \int_0^\infty \mathrm{Ai}(z+u) \frac{du}{\sqrt{u}} = 2^{2/3} \pi \mathrm{Ai}^2\left(\frac{z}{2^{2/3}}\right) \equiv I(z) \quad (210) $$

and differentiate it twice with respect to $z$. Using the Airy differential equation $\mathrm{Ai}''(z) = z\mathrm{Ai}(z)$, one obtains

$$ \begin{aligned} \int_0^\infty du \sqrt{u} \mathrm{Ai}(z+u) &= I''(z) - zI(z) \\ &= \pi 2^{1/3} \left\{ \left[ \mathrm{Ai}'\left(\frac{z}{2^{2/3}}\right) \right]^2 - \frac{z}{2^{2/3}} \mathrm{Ai}^2\left(\frac{z}{2^{2/3}}\right) \right\}. \end{aligned} \quad (211) $$

It then follows from Eq. (207), upon setting $d = 1$, that

$$ F_1(z) = \mathrm{Ai}'^2(z) - z\mathrm{Ai}^2(z), \quad (212) $$

thus recovering the result obtained in Sec. III A [see Eq. (31)], which coincides with the well-known RMT result [49,50]. One obtains similar quadratic forms in $\mathrm{Ai}(z)$ and $\mathrm{Ai}'(z)$ with polynomial coefficients in $z$ in any odd space dimension by repeated application of the Airy operator $(\partial_z^2 - z)$ on $I(z)$. For instance, in $d = 3$

$$ F_3(z) = \frac{1}{12\pi} [2z^2 \mathrm{Ai}(z)^2 - \mathrm{Ai}(z)\mathrm{Ai}'(z) - 2z\mathrm{Ai}'(z)^2]. \quad (213) $$

In $d = 2$ one can use the Airy equation and find

$$ F_2(z) = \frac{1}{2^{8/3}\pi} \left[ -\mathrm{Ai}'(2^{1/3}z) - 2^{1/3}z\mathrm{Ai}_1(2^{1/3}z) \right], \quad (214) $$

where $\mathrm{Ai}_1(z) = \int_z^\infty dx \mathrm{Ai}(x)$.

Asymptotic behavior of $F_d(z)$. Here, we give the asymptotic behavior of the scaling functions $F_d(z)$, the full form of which are plotted in Fig. 8 for $d = 1, 2, 3$. We first consider the $z \to +\infty$ limit. In this limit, the Airy function has the leading asymptotic behavior [80]

$$ \mathrm{Ai}(z) \sim \frac{1}{2\sqrt{\pi}} z^{-1/4} \exp\left(-\frac{2}{3}z^{3/2}\right). \quad (215) $$

Substituting this asymptotic behavior in Eq. (207), expanding for large $z$, one gets to leading order

$$ F_d(z) \approx (8\pi)^{-d/2} z^{-d/4} \exp\left(-\frac{4}{3} z^{3/2}\right) \quad \text{as } z \to \infty. \quad (216) $$

For the other side $z \to -\infty$, it is more convenient to use the representation in Eq. (203). We set $z = -|z|$ and scale $\tau |z| = t$. This makes the order $\tau^3$ term to be $|z|^{-3} t^3/12$ which can then be dropped for large $|z|$. The resulting Bromwich

<figure><img src="image_9.png" /><figcaption>FIG. 8. Plot of the scaling functions $F_d(z)$ in Eq. (207) 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>

contour integral can be easily evaluated to give the leading asymptotic behavior

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

One can show that when $z \to -\infty$, i.e., when $r \ll r_{\mathrm{edge}}$, the asymptotic behavior in Eq. (217) matches smoothly with the bulk density given in Eq. (183). This can be seen by noting that one can write $r_{\mathrm{edge}} = 1/(2w_N^3\alpha^4)$ and then writing $r = r_{\mathrm{edge}} - |z|w_N$ in Eq. (183), which then becomes

$$ w_N^d \rho(r_{\mathrm{edge}} - |z|w_N) \approx \frac{|z|^{d/2}}{(4\pi)^{d/2} \Gamma(1 + d/2)}, \quad (218) $$

which coincides with the behavior in Eq. (217).

## 2. Kernel near the edge

In order to analyze the full kernel $K_\mu(\mathbf{x}, \mathbf{y})$ near the edge, we introduce the following scaled dimensionless coordinates near a point $\mathbf{r}_{\mathrm{edge}}$ on the circle of radius $r_{\mathrm{edge}}$ as

$$ \mathbf{x} = \mathbf{r}_{\mathrm{edge}} + w_N \mathbf{a}, \quad \mathbf{y} = \mathbf{r}_{\mathrm{edge}} + w_N \mathbf{b}, \quad (219) $$

where $w_N$ is given in Eq. (202) and

$$ r_{\mathrm{edge}} = \sqrt{2\mu/(m\omega^2)} \approx [\Gamma(d+1)]^{1/(2d)} N^{1/(2d)} \sqrt{2}/\alpha. \quad (220) $$

Following the analysis for the edge density in the previous section, we insert (219) into the expansion of $K_\mu(\mathbf{x}, \mathbf{y})$ up to order $O(t^3)$ given in Eq. (174). We note that if one takes $|\mathbf{a}|, |\mathbf{b}| = O(1)$, the diffusion part in Eq. (174) scales as

$$ \frac{m}{2\pi\hbar t} (\mathbf{x} - \mathbf{y})^2 \sim \frac{w_N^2}{t} \sim \frac{N^{-1/(3d)}}{t} = O(1), \quad (221) $$

where we have used Eq. (199). And the analysis for the other terms is similar to the one performed for the density since one has

$$ \mathbf{x}^2 + \mathbf{y}^2 = 2 r_{\mathrm{edge}}^2 [1 + w_N(a_n + b_n)] + O(w_N^2), \quad (222) $$

$$ (\mathbf{x} - \mathbf{y})^2 = w_N^2 (\mathbf{a} - \mathbf{b})^2 \ll \mathbf{x}^2 + \mathbf{y}^2, \quad (223) $$
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 the projections of $\mathbf{a}$ and $\mathbf{b}$ in the radial direction. Putting all together gives

$$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^{3/2}\tau} - \frac{(a_n+b_n)\tau}{2^{1/3}} + \frac{\tau^3}{3}}, \quad (224)$$

with $C_d = (2^{4/3}\sqrt{\pi})^d$ and where we have rescaled $t = N^{-1/(3d)}/(\omega[\Gamma(d+1)]^{1/(3d)})$ $\tau$, followed by $\tau \to 2^{2/3}\tau$ as in the previous section.

One can make a further simplification of Eq. (224) by using the integral representation of the diffusive propagator

$$\frac{e^{-\frac{(\mathbf{a}-\mathbf{b})^2}{4D\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 (225)$$

We choose $D = 2^{2/3}$ and use this in Eq. (224). This gives our final result for the scaling behavior of the edge kernel,

$$K_{\mu}(\mathbf{x}, \mathbf{y}) \approx \frac{1}{w_N^d} \mathcal{K}_d^{\text{edge}} \left( \frac{\mathbf{x} - \mathbf{r}_{\text{edge}}}{w_N}, \frac{\mathbf{y} - \mathbf{r}_{\text{edge}}}{w_N} \right), \quad (226)$$

where the scaling function is given explicitly by

$$\mathcal{K}_d^{\text{edge}}(\mathbf{a}, \mathbf{b}) = \int \frac{d^d q}{(2\pi)^d} e^{-i \mathbf{q} \cdot (\mathbf{a} - \mathbf{b})} \text{Ai}_1 \left( 2^{\frac{2}{3}} q^2 + \frac{a_n + b_n}{2^{1/3}} \right), \quad (227)$$

while the function $\text{Ai}_1(z)$ is given by

$$\text{Ai}_1(z) = \int_{\Gamma} \frac{d\tau}{2\pi i} \frac{1}{\tau} e^{-z\tau + \tau^{3/3}} = \int_{z}^{\infty} \text{Ai}(u) du. \quad (228)$$

This edge kernel is a novel result and generalizes the standard Airy kernel in $d=1$ to higher dimensions. Indeed, putting $d=1$ in Eq. (227) we get

$$\mathcal{K}^{\text{edge}}(a,b) = \int_{-\infty}^{\infty} \frac{dq}{2\pi} e^{iq(a-b)} \int_{2^{2/3}q^2 + 2^{-1/3}(a+b)}^{\infty} \text{Ai}(z) dz. \quad (229)$$

Making a shift $z = 2^{2/3}q^2 + 2^{-1/3}(a+b) + u$ gives

$$\mathcal{K}^{\text{edge}}(a,b) = \int_{-\infty}^{\infty} \frac{dq}{2\pi} e^{iq(a-b)} \times \int_{0}^{\infty} \text{Ai}[u + 2^{2/3}q^2 + 2^{-1/3}(a+b)] du. \quad (230)$$

We next use a nontrivial identity involving Airy functions [80]

$$\int_{-\infty}^{\infty} \frac{dq}{2\pi} e^{-iq(v-v')} \text{Ai}[2^{2/3}q^2 + 2^{-1/3}(v+v')] \\ = 2^{-\frac{2}{3}} \text{Ai}(v)\text{Ai}(v'). \quad (231)$$

Choosing $v = a + 2^{-2/3}u$ and $v' = b + 2^{-2/3}u$, substituting this identity in Eq. (230) and rescaling $u \to 2^{-2/3}u$ gives

$$\mathcal{K}^{\text{edge}}(a,b) = \int_{0}^{\infty} du \, \text{Ai}(a+u) \text{Ai}(b+u). \quad (232)$$

Since $\text{Ai}(z)$ satisfies the differential equation $\text{Ai}''(z) - z\text{Ai}(z) = 0$ we replace $\text{Ai}(z)$ by $\text{Ai}''(z)/z$ in Eq. (232). Next,

we use the identity

$$\frac{1}{(u+a)(u+b)} = \frac{1}{b-a} \left[ \frac{1}{u+a} - \frac{1}{u+b} \right] \quad (233)$$

and integrate by parts. This then reduces Eq. (232) to the standard Airy kernel form

$$\begin{aligned} \mathcal{K}^{\text{edge}}(a,b) &= K_{\text{Airy}}(a,b) \\ &= [\text{Ai}(a) \text{Ai}'(b) - \text{Ai}'(a) \text{Ai}(b)]/(a-b). \end{aligned} \quad (234)$$

## VII. GENERAL D-DIMENSIONAL SOFT POTENTIALS AT FINITE T AND UNIVERSALITY

In this section, we finally consider the most general case of arbitrary $d$ and finite $T$. In addition, we show that the results obtained for the harmonic oscillator in the previous sections can be extended to a very broad class of smooth confining potentials. We start by generalizing the method based on the Euclidean propagator, introduced in the previous section, to finite temperature (Sec. VII A). To study the large-$N$ limit, we only need the small time expansion of the propagator, which is obtained for a general potential in Sec. VII B. Using this expansion, we obtain results in the bulk in Sec. VII C and at the edge of the Fermi gas in Sec. VII D. Most of our formulas will be valid for arbitrary confining smooth potentials $V(\mathbf{x})$. For the spherically symmetric potentials, in particular those of the type

$$V(\mathbf{x}) = V(r) = V_0 \left(\frac{r}{r_0}\right)^p, \quad (235)$$

we will obtain additional explicit results. Note that this parametrization contains some arbitrariness but final results only depend on $V_0/r_0^p$. For convenience, we choose

$$V_0 = \frac{\hbar^2}{2mr_0^2} \quad (236)$$

so that $V_0, r_0$ are, respectively, of the order of the single-particle ground-state energy and its radius. The harmonic oscillator is recovered for $p=2, r_0 = 1/\alpha$, and $V_0 = \frac{1}{2}\hbar\omega$.

### A. Representation of the finite-T kernel using the Euclidean propagator

To deal with finite $T$ we need to generalize to arbitrary $d$ the method explained in detail for $d=1$ in Sec. IV. There we showed that the canonical and grand-canonical ensembles lead to the same results for the $n$-point correlations for large $N$. This allows us to work here in the grand-canonical ensemble, where the method of the one-body Euclidean propagator can also be applied.

To study the problem at finite $T$, one considers arbitrary excited eigenstates of the $N$-body Hamiltonian $\hat{H}_N$, with arbitrary potential $V(\mathbf{x})$ in Eq. (1). These states are labeled by a set of occupation numbers $n_k = 0, 1$. As discussed in Sec. IV, to each excited state one associates a kernel given by

$$K(\mathbf{x}, \mathbf{y}; \{n_k\}) = \sum_{\mathbf{k}} n_{\mathbf{k}} \psi_{\mathbf{k}}^{*}(\mathbf{x}) \psi_{\mathbf{k}}(\mathbf{y}), \quad (237)$$

which generalizes the $d=1$ formula (59). We recall that the $\epsilon_{\mathbf{k}}$ and $\psi_{\mathbf{k}}$ are the eigenenergies and eigenfunctions of
the one-body Hamiltonian. Following the steps presented in Sec. IV we obtain the generalization of the set of formulas (91)–(95). In particular, the n-point correlation functions $R_n$ are given, in any d, by determinants constructed from the grand-canonical kernel, as in formula (93). This kernel is obtained by performing the average over the $n_k$ in the grand-canonical ensemble, leading to

$$ K_{\tilde{\mu}}(\mathbf{x}, \mathbf{y}) := \langle K(\mathbf{x}, \mathbf{y}; \{n_k\}) \rangle \\ = \sum_{\mathbf{k}} \frac{1}{1 + \exp[\beta(\epsilon_{\mathbf{k}} - \tilde{\mu})]} \psi_{\mathbf{k}}^{*}(\mathbf{x})\psi_{\mathbf{k}}(\mathbf{y}), \quad (238) $$

where the chemical potential $\tilde{\mu}$ is related to N and the temperature T via the relation

$$ N = \sum_{\mathbf{k}} \frac{1}{1 + \exp[\beta(\epsilon_{\mathbf{k}} - \tilde{\mu})]}. \quad (239) $$

We recall that, as in $d=1$, we reserve the notations $\mu$ and $K_\mu(\mathbf{x}, \mathbf{y})$ for the zero-temperature chemical potential and kernel, respectively, while we denote $\tilde{\mu}$ and $K_{\tilde{\mu}}(\mathbf{x}, \mathbf{y})$ their finite-temperature counterparts. We recall that, as the temperature goes to zero, $\tilde{\mu} \to \mu$, but at finite temperature, not only the chemical potential $\tilde{\mu}$ differs from $\mu$, but also the full kernel functions are different.

A useful representation of the finite-temperature kernel can now be derived in terms of the Euclidean propagator associated to the one-body Hamiltonian (1) with arbitrary $V(\mathbf{x})$, defined in Eq. (166). Using the formula (164) for the zero-temperature kernel $K_\mu(\mathbf{x}, \mathbf{y})$, one can rewrite

$$ K_{\tilde{\mu}}(\mathbf{x}, \mathbf{y}) = \int_{0}^{\infty} d\mu \frac{\partial K_{\mu}(\mathbf{x}, \mathbf{y})}{\partial \mu} \frac{1}{1 + \exp[\beta(\mu - \tilde{\mu})]}. \quad (240) $$

Now, using Eq. (168) we obtain the following representation of the finite-temperature kernel:

$$ K_{\tilde{\mu}}(\mathbf{x}, \mathbf{y}) = \int_{0}^{\infty} d\mu' \frac{1}{1 + \exp[\beta(\mu' - \tilde{\mu})]} \times \int_{\Gamma} \frac{dt}{2\pi\hbar i} \exp\left(\frac{\mu't}{\hbar}\right) G(\mathbf{x}, \mathbf{y}; t), \quad (241) $$

where we have used $\mu'$ as a running integration variable, rather than $\mu$, which in the remaining will always denote the chemical potential at zero temperature.

## B. Small time expansion of the propagator for generic potential

We have seen in the previous section, for the d-dimensional simple harmonic oscillator at $T=0$, that in the limit of large $\mu$ the kernel and density can be extracted from the short-time behavior of the Euclidean propagator. As shown below, this remains true even at finite temperature, in the relevant regimes studied here. This extends the short-time expansion analysis initiated in Ref. [81]. However, for the sake of completeness, we now display the result for a general soft potential $V(\mathbf{x})$, and provide an independent probabilistic derivation of this expansion in Appendix A.

The expansion of $G(\mathbf{x}, \mathbf{y}; t)$ for a soft potential is, to $O(t^3)$, given by

$$ G(\mathbf{x}, \mathbf{y}; t) \sim \left(\frac{m}{2\pi\hbar t}\right)^{\frac{d}{2}} \exp\left[-\frac{m}{2\hbar t}(\mathbf{x} - \mathbf{y})^2\right] \times \exp\left[-\frac{t}{\hbar}S_1(\mathbf{x}, \mathbf{y}) - \frac{t^2}{2m}S_2(\mathbf{x}, \mathbf{y}) + \frac{t^3}{2m\hbar}S_3(\mathbf{x}, \mathbf{y})\right], \quad (242) $$

where

$$ S_1(\mathbf{x}, \mathbf{y}) = \int_0^1 du V[\mathbf{x} + u(\mathbf{y} - \mathbf{x})], \quad (243) $$

$$ S_2(\mathbf{x}, \mathbf{y}) = \int_0^1 du u(1-u) \nabla^2 V[\mathbf{x} + u(\mathbf{y} - \mathbf{x})], \quad (244) $$

$$ S_3(\mathbf{x}, \mathbf{y}) = \int_0^1 du \int_0^1 dw [\min(u, w) - uw] \times \nabla V[\mathbf{x} + u(\mathbf{y} - \mathbf{x})] \cdot \nabla V[\mathbf{x} + w(\mathbf{y} - \mathbf{x})] - \frac{\hbar^2}{4m} \int_0^1 du u^2 (1-u)^2 \nabla^2 \nabla^2 V[\mathbf{x} + u(\mathbf{y} - \mathbf{x})]. \quad (245) $$

With some work it can be checked that the above result agrees with the short-time expansion given in Eq. (172) when applied to the simple harmonic oscillator. We see that the occurrence of derivatives of first and higher order in $V$, resulting from the expansion of the term $V[\mathbf{x}(1-u) + \mathbf{y}u + \sqrt{D_0t}\mathbf{B}_u]$ in Eq. (A1), where $\mathbf{B}_u$ denotes the d-dimensional Brownian bridge, requires that the potential is sufficiently smooth within the neighborhood of the direct path between $\mathbf{x}$ and $\mathbf{y}$. At coinciding points $\mathbf{x} = \mathbf{y}$ (the case where one computes the density), one finds

$$ S_1(\mathbf{x}, \mathbf{x}) = V(\mathbf{x}), \quad (246) $$

$$ S_2(\mathbf{x}, \mathbf{x}) = \frac{1}{6} \nabla^2 V(\mathbf{x}), \quad (247) $$

$$ S_3(\mathbf{x}, \mathbf{x}) = \frac{1}{12} [\nabla V(\mathbf{x})]^2 - \frac{\hbar^2}{120m} \nabla^2 \nabla^2 V(\mathbf{x}). \quad (248) $$

Using these results, we now analyze successively the bulk and the edge regimes.

### C. Results in the bulk

#### 1. Zero temperature

The calculation of the bulk density and kernel follows exactly those of the simple harmonic oscillator in Secs. VIB 1 and VIB 2. Here, only the term to $O(t)$ is retained and one obtains

$$ K_{\mu}(\mathbf{x}, \mathbf{y}) = \theta[\mu - S_1(\mathbf{x}, \mathbf{y})] \left\{ \frac{[\mu - S_1(\mathbf{x}, \mathbf{y})]m}{2\pi^2 \hbar^2 (\mathbf{x} - \mathbf{y})^2} \right\}^{\frac{d}{2}} \times J_{\frac{d}{2}} \left\{ \sqrt{\frac{2m(\mathbf{x} - \mathbf{y})^2[\mu - S_1(\mathbf{x}, \mathbf{y})]}{\hbar^2}} \right\}. \quad (249) $$

Let us first discuss the normalized particle density is given by $\rho_N(\mathbf{x}) = K_\mu(\mathbf{x}, \mathbf{x})/N$. Again, we see that its support
$S_1(\mathbf{x}, \mathbf{y}) = V(\mathbf{x}) < \mu$ is finite and using Eq. (246) is obtained as

$$ \rho_N(\mathbf{x}) = \left(\frac{m}{2\hbar^2\pi}\right)^{\frac{d}{2}} \frac{\theta[\mu - V(\mathbf{x})]}{N\Gamma(1 + \frac{d}{2})} [\mu - V(\mathbf{x})]^{\frac{d}{2}}. \quad (250) $$

From this formula, by integrating over $\mathbf{x}$, one can calculate the Fermi energy (or temperature) $T_F = \mu$, from the condition $\int d\mathbf{x} \rho_N(\mathbf{x}) = 1$.

For explicit calculations, let us focus on the case of spherically symmetric potential of the form (235) with $p > 0$. The density vanishes at the edge at $r = r_{\text{edge}}$ such that

$$ V(r_{\text{edge}}) = \mu \iff r_{\text{edge}} = r_0 \left(\frac{\mu}{V_0}\right)^{1/p}. \quad (251) $$

For the harmonic oscillator in $d=1$ one recovers $r_{\text{edge}} = \sqrt{2N}/\alpha$ using that in that case $\mu = N\hbar\omega$. In the case of general $p$, we obtain

$$ T_F = \mu = a_{p,d} V_0 N^{\frac{2p}{d(p+2)}} \quad (252) $$

with $a_{p,d} = (4\pi)^{\frac{p}{2+p}} \{p\Gamma(1+d/2)/[S_d B(1+d/2, d/p)]\}^{2p/(d(p+2))}$ where $S_d = 2\pi^{d/2}/\Gamma(d/2)$ and $B(p,q) = \int_0^1 u^{p-1}(1-u)^{q-1} du$ and where we used (236). Consequently, one obtains that

$$ r_{\text{edge}} \sim N^{2/(d(p+2))} \quad (253) $$

for large $N$.

Consider now the kernel given in Eq. (249). As in the case of the simple harmonic oscillator, we consider two generic points $\mathbf{x}$ and $\mathbf{y}$ in the bulk (far from the edges) close to each other, with a separation of order $|\mathbf{x} - \mathbf{y}| \sim [N\rho_N(\mathbf{x})]^{-1/d}$. Equation (249) then simplifies to the scaling form

$$ K_{\mu}(\mathbf{x}, \mathbf{y}) \approx \ell^{-d}(\mathbf{x}) K_{\mu}^{\text{bulk}}[|\mathbf{x} - \mathbf{y}|/\ell(\mathbf{x})], \quad (254) $$

where $\ell(\mathbf{x}) = [N\rho_N(\mathbf{x})\gamma_d]^{-1/d}$ is the typical separation between fermions in the bulk, $\gamma_d = \pi^{d/2}[\Gamma(d/2 + 1)]$, and $K_{\mu}^{\text{bulk}}$ is the same scaling function, as given in Eq. (194), respectively, for the harmonic oscillator. The dependence on the potential $V(\mathbf{x})$ thus enters only through the local density $\rho_N(\mathbf{x})$ via the scale factor $\ell(\mathbf{x})$. However, the scaling function associated with the bulk kernel in Eq. (194) is completely universal for all $V(\mathbf{x})$.

The above result is the general form of the LDA [2], which is normally obtained from semiclassical or physical arguments. The range of validity of this approximation can, as was the case for the simple harmonic oscillator, be established by examining the corrections due to the quadratic and cubic terms in $t$ in the short-time expansion of the propagator. Here, we define the two-point energy function

$$ \epsilon_2(\mathbf{x}, \mathbf{y}) = \mu - S_1(\mathbf{x}, \mathbf{y}), \quad (255) $$

the range of validity of the LDA approximation is thus given by

$$ \epsilon_2(\mathbf{x}, \mathbf{y}) \gg \left[ \hbar^2 \frac{S_2(\mathbf{x}, \mathbf{y})}{m} \right]^{\frac{1}{2}} \quad (256) $$

and

$$ \epsilon_2(\mathbf{x}, \mathbf{y}) \gg \left[ \hbar^2 \frac{S_3(\mathbf{x}, \mathbf{y})}{m} \right]^{\frac{1}{3}}. \quad (257) $$

If we consider a trapped system with a potential which grows in a polynomial fashion we see, from power counting, that for large values of $\mathbf{x}$ and $\mathbf{y}$, the second of the above inequalities determines the validity of the LDA.

## 2. Bulk statistics at finite temperature and $T \sim T_F$

At zero temperature, there is a unique length scale associated to the quantum fluctuation in the confining potential, denoted by $r_0$ for the class of potentials (235) above. By contrast, at finite temperature, for a generic confining potential, there are two natural length scales in the problem. The first one is the thermal de Broglie wavelength

$$ \lambda_T = \hbar \sqrt{\frac{2\pi}{mT}} \quad (258) $$

of the free fermions. The second length is set by the temperature and the confining potential and is purely classical. In the case of a power-law spherically symmetric potential (235) it is given by

$$ \beta V(r_T) = 1 \iff r_T = r_0 \left(\frac{T}{V_0}\right)^{1/p}. \quad (259) $$

Now, let us discuss the energy scales involved. At zero temperature the natural energy scale is the Fermi energy $\mu = T_F$ defined and discussed in the previous section, where we showed that $\mu$ is large for large $N$, i.e., $\mu \sim N^{2p/(d(p+2))}$ for the class of potentials defined in Eq. (235). Significant changes from the $T=0$ properties are expected in the bulk only when $T \sim T_F$. In what follows, we will focus on this bulk regime and consider $T \sim T_F = \mu$. Consequently, the fugacity which we denote as

$$ \zeta = \exp(\beta \tilde{\mu}) \quad (260) $$

is also of $O(1)$ when $N$ is large. Introducing the variable $u = \beta\mu'$ in Eq. (241) and making the change of variable $t = \tau\beta\hbar$ we find

$$ K_{\tilde{\mu}}(\mathbf{x}, \mathbf{y}) = \int_{0}^{\infty} du \frac{1}{1 + \zeta^{-1} \exp(u)} \times \int_{\Gamma} \frac{d\tau}{2\pi i} \exp(u\tau) G(\mathbf{x}, \mathbf{y}; \beta\hbar\tau), \quad (261) $$

which at this stage is still an exact equation.

We now analyze this equation in the large-$N$ limit and the regime $T \sim T_F$. This means that $\beta$ is small, hence the time variable $t$ in the above propagator $G(\mathbf{x}, \mathbf{y}; t = \beta\hbar\tau)$ is also small. Hence, we can use the short-time expansion of the Euclidean propagator, as was done for the harmonic oscillator case. As discussed, there to describe the bulk, one only needs to keep the first term of order $O(t)$ and one obtains

$$ K_{\tilde{\mu}}(\mathbf{x}, \mathbf{y}) = \int_{0}^{\infty} du \frac{1}{1 + \zeta^{-1} \exp(u)} \int_{\Gamma} \frac{d\tau}{2\pi i} \left( \frac{m}{2\pi \tau \beta \hbar^2} \right)^{\frac{d}{2}} \times \exp \left[ - \frac{\pi (\mathbf{x} - \mathbf{y})^2}{\tau \lambda_T^2} + \tau u - \tau \beta S_1(\mathbf{x}, \mathbf{y}) \right], \quad (262) $$

where $\lambda_T$ is given in Eq. (258), $\zeta$ by (260) where $\tilde{\mu}$ is fixed as a function of $N$ by Eq. (239).
To obtain the density, we start by considering the kernel at coinciding points

$$K_{\bar{\mu}}(\mathbf{x}, \mathbf{x}) = \int_{0}^{\infty} du \frac{1}{1 + \zeta^{-1}\exp(u)} \int_{\Gamma} \frac{d\tau}{2\pi i} \left(\frac{m}{2\pi \tau \beta \hbar^2}\right)^{\frac{d}{2}} \exp[\tau u - \tau \beta V(\mathbf{x})]. \quad (263)$$

The appropriate bulk scaling is thus to choose $\mathbf{x}$ such that $\beta V(\mathbf{x}) \sim O(1)$. Performing the explicit integral over the Bromwich contour using (181) we obtain

$$K_{\bar{\mu}}(\mathbf{x}, \mathbf{x}) = \left(\frac{m}{2\pi \beta \hbar^2}\right)^{\frac{d}{2}} \int_{0}^{\infty} du \frac{1}{1 + \zeta^{-1}\exp(u)} \frac{\theta[u - \beta V(\mathbf{x})][u - \beta V(\mathbf{x})]^{\frac{d}{2}-1}}{\Gamma(\frac{d}{2})}. \quad (264)$$

This can then be rewritten, shifting the integral over $u$, and leads to our main result for the density

$$\rho_N(\mathbf{x}) = \frac{1}{N} K_{\bar{\mu}}(\mathbf{x}, \mathbf{x}) = \frac{1}{N\Gamma(\frac{d}{2})\lambda_T^d} \int_0^\infty dq \frac{q^{\frac{d}{2}-1}}{1 + \zeta^{-1}\exp[q + \beta V(\mathbf{x})]} = -\frac{1}{N\lambda_T^d} \mathrm{Li}_{d/2}(-\zeta e^{-\beta V(\mathbf{x})}), \quad (265)$$

where $\lambda_T$ is the thermal wavelength of the free fermions given in Eq. (258) and $\mathrm{Li}_n(x) = \sum_{k=1}^{\infty} x^k/k^n$ is the polylogarithm function.

Let us provide an explicit example in the case of the power-law potentials (235). We define the scaling variables, which generalize the 1d result given in Eqs. (99) and (100)

$$y = \left(\frac{T_F}{T}\right)^d, \quad z = r/r_T, \quad (266)$$

where $T_F$ is given in Eq. (252) and $r_T$ is the classical thermal confining length introduced in Eq. (259). The explicit calculation predicts that the density is given by eliminating the fugacity $\zeta$ in the following pair of equations:

$$\rho_N(\mathbf{x}) = -\frac{B_2}{\sqrt{yN^{\frac{2}{p+2}}}} \mathrm{Li}_2(-\zeta e^{-z^p}), \quad (267)$$

$$1 = -\frac{B_1}{y^{\frac{p+2}{2p}}} \mathrm{Li}_{\frac{d}{2}+\frac{d}{p}}(-\zeta), \quad (268)$$

where $B_2 = (4\pi)^{-d/2}r_0^{-d}a_{p,d}^{d/2}$, $a_{p,d}$ is given below Eq. (252), and $B_1 = \Gamma(1 + d/2 + d/p)$. Using the identity

$$\int_{0}^{+\infty} dr \, r^a \mathrm{Li}_b(-\zeta e^{-r^p}) = \frac{1}{p} \Gamma\left(\frac{a+1}{p}\right) \mathrm{Li}_{\frac{1+a+b}{p}}(-\zeta), \quad (269)$$

we can check that the Eq. (268) is the normalization condition $\int d^d x \rho_N(\mathbf{x}) = 1$ of the density given by Eq. (267). Setting $p=2$ in the above formula gives the result for the $d$-dimensional harmonic oscillator. In particular in $d=1$, using that $\mathrm{Li}_1(-\zeta) = -\ln(1+\zeta)$, the implicit equation can be solved and one recovers the bulk scaling function for the density presented in Eqs. (109) and (110) in Sec. V.

The above result agrees with the LDA approximation which we briefly recall here for completeness. In the LDA, the local (unnormalized density) in position and momentum is given by the Fermi factor

$$n(\mathbf{x}, \mathbf{p}) = \frac{1}{1 + \exp[\beta E(\mathbf{x}, \mathbf{p}) - \beta \bar{\mu}]}, \quad (270)$$

where $E(\mathbf{x}, \mathbf{p})$ is given by the classical energy

$$E(\mathbf{x}, \mathbf{p}) = \frac{\mathbf{p}^2}{2m} + V(\mathbf{x}). \quad (271)$$

Integration over the momentum with the appropriate Planck volume normalization then gives the local density in space as

$$N \rho_N(\mathbf{x}) = K_{\bar{\mu}}(\mathbf{x}, \mathbf{x}) = \frac{1}{h^d} \int \frac{d\mathbf{p}}{1 + \exp[\beta E(\mathbf{p}, \mathbf{x}) - \beta \bar{\mu}]}. \quad (272)$$

Carrying out the angular integration in Eq. (272) yields

$$\rho_N(\mathbf{x}) = \frac{2\pi^{\frac{d}{2}}}{N\Gamma(\frac{d}{2})h^d} \int_0^\infty \frac{dp \, p^{d-1}}{1 + \exp[\beta E(p, \mathbf{x}) - \beta \bar{\mu}]}. \quad (273)$$

Finally, making the substitution $q = \beta p^2/2m$ shows that the LDA approximation (273) agrees with Eq. (265).

The formula for the bulk kernel at two unequal points is obtained by analyzing Eq. (262) in the regime of separation $|\mathbf{x} - \mathbf{y}| \sim \lambda_T$, and of temperature such that $\beta V(\mathbf{x}) = O(1)$. In this limit, by expanding $V$ in the integrand of (243), we can approximate $S_1(\mathbf{x}, \mathbf{y})$ by its leading term $V(\mathbf{x})$, using that $\beta \lambda_T |\nabla V(\mathbf{x})| \ll 1$ (always valid inside the bulk in this temperature
regime, as can be checked explicitly). Replacing $S_1(\mathbf{x}, \mathbf{y})$ by $V(\mathbf{x})$ in Eq. (262) and using the change of variable $u - \beta V(\mathbf{x}) = q$ and the integral representation of the Bessel function in Eq. (190), we arrive at

$$ K_{\tilde{\mu}}(\mathbf{x}, \mathbf{y}) = \frac{1}{\pi^{\frac{d-2}{4}} \lambda_T^d} \left( \frac{\lambda_T}{|\mathbf{x} - \mathbf{y}|} \right)^{\frac{d}{2}-1} \int_0^{\infty} dq \frac{q^{\frac{d-2}{4}}}{1 + \zeta^{-1} \exp[+\beta V(\mathbf{x})]} J_{\frac{d}{2}-1} \left( 2\sqrt{q\pi} \frac{|\mathbf{x} - \mathbf{y}|}{\lambda_T} \right). \quad (274) $$

This formula crosses over to the zero-temperature kernel (194) when $\lambda_T \gg \ell(\mathbf{x})$ where $\ell(\mathbf{x})$ is the typical separation between particles defined in Eq. (193). This can be seen by performing a change of variable $q \to \beta q$, such that the Fermi factor becomes a Heaviside step function. Thus, in practice the above formula is valid in the range of separation $|\mathbf{x} - \mathbf{y}| \sim \min(\lambda_T, \ell(\mathbf{x}))$. Beyond this scale the kernel decays to zero. Note that in $d=1$ and $3$ the formula simplifies slightly since one has, respectively, $J_{-\frac{1}{2}}(x) = \sqrt{2/\pi x} \cos(x)$ and $J_{\frac{1}{2}}(x) = \sqrt{2/\pi x} \sin(x)$.

### D. Results at the edge

In the edge region we expect the effects of fluctuations to be larger than in the bulk. To study the bulk we had to scale the temperature as $T \sim T_F$. By contrast, in this section, the relevant regime will involve lower temperature, $T \ll T_F$, hence, everywhere the variable $\beta\mu$ will be considered to be large. To estimate the chemical potential $\tilde{\mu}$ in this range of temperature, we can thus use the Sommerfeld expansion valid for large $\beta\mu \gg 1$ [3]:

$$ \mu - \tilde{\mu} = \frac{\pi^2}{6} \frac{1}{\beta^2} \frac{\tilde{\rho}'(\mu)}{\tilde{\rho}(\mu)} + O\left[ \frac{1}{(\beta\mu)^4} \right], \quad (275) $$

where we denote by $\tilde{\rho}(\epsilon) = \sum_k \delta(\epsilon - \epsilon_k)$ the density of states (in energy). Note that for the harmonic oscillator in $d=1$ the density of states is constant, and the corrections are exponentially small.

In this section, we consider the kernel for a pair of points $\mathbf{x}, \mathbf{y}$ both located near a point $\mathbf{r}_{\text{edge}}$ at the edge, at finite temperature, in arbitrary $d$ and for the large class of confining potential $V(r)$ described at the beginning of Sec. VII. To proceed, we set

$$ \mathbf{x} = \mathbf{r}_{\text{edge}} + \mathbf{a}', \quad \mathbf{y} = \mathbf{r}_{\text{edge}} + \mathbf{b}', \quad (276) $$

where $|\mathbf{r}_{\text{edge}}| = r_{\text{edge}}$ is defined in Eq. (251) and we assume $|\mathbf{a}'|, |\mathbf{b}'| \ll r_{\text{edge}}$. The goal of this section is to show that the properly centered and scaled edge kernel becomes universal, i.e., independent of the details of the potential.

Substituting (276) in Eqs. (243)–(245) and expanding the terms $S_1, S_2$, and $S_3$, in a gradient expansion, we obtain

$$ S_1 = V(\mathbf{r}_{\text{edge}}) + \frac{1}{2} \nabla V(\mathbf{r}_{\text{edge}}) \cdot (\mathbf{a}' + \mathbf{b}') + \dots, \quad (277) $$

$$ S_2 = \frac{1}{6} \nabla^2 V(\mathbf{r}_{\text{edge}}) + \frac{1}{12} \nabla [\nabla^2 V(\mathbf{r}_{\text{edge}})] \cdot (\mathbf{a}' + \mathbf{b}') + \dots, \quad (278) $$

$$ S_3 = \frac{1}{12} [\nabla V(\mathbf{r}_{\text{edge}})]^2 - \frac{\hbar^2}{120m} \nabla^2 \nabla^2 V(\mathbf{r}_{\text{edge}}) + \dots, \quad (279) $$

Note that for potentials increasing as $V(r) \sim r^p$ for large $r$, each derivative brings an additional factor $1/r_{\text{edge}}$ in the expansion, hence, it is not necessary to keep higher-order terms, at least for the potentials of this class.

Let us start with Eq. (241) and substitute in it the short-time expansion (245). We will justify a posteriori under what conditions the short-time expansion can be stopped at order $O(t^3)$. We make a change of variable $\beta(\mu' - \tilde{\mu}) = -bu$ in Eq. (241), where for convenience we have introduced a dimensionless parameter $b$ whose value will be chosen later. Setting $\mathbf{x}$ and $\mathbf{y}$ both close to $\mathbf{r}_{\text{edge}}$, as in Eq. (276), and using (277)–(279) we obtain

$$ K_{\tilde{\mu}}(\mathbf{x}, \mathbf{y}) \approx \frac{b}{\beta} \int_{-\infty}^{\infty} \frac{du}{1 + \exp(-bu)} \int_{\Gamma} \frac{dt}{2\pi\hbar i} \left( \frac{m}{2\pi\hbar t} \right)^{\frac{d}{2}} \times \exp \left[ -\frac{m}{2\hbar t} (\mathbf{a}' - \mathbf{b}')^2 \right] \exp \left[ -\frac{ubt}{\beta\hbar} - \frac{t}{2\hbar} |\nabla V(\mathbf{r}_{\text{edge}})| (a'_n + b'_n) + \frac{t^3}{24m\hbar} |\nabla V(\mathbf{r}_{\text{edge}})|^2 \right], \quad (280) $$

where the upper bound $\tilde{\mu}\beta/b$ has been replaced by $+\infty$ since we are studying the limit of large $\tilde{\mu}$ (large $N$). In deriving this equation we have kept both terms of $S_1$, but neglected $S_2$ and the second term in $S_3$ in Eqs. (277)–(279), which will be justified later. Here, $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 as in Sec. VI C. The term linear in time $t/\hbar[\tilde{\mu} - V(\mathbf{r}_{\text{edge}})] = t/\hbar(\tilde{\mu} - \mu) \sim \hbar\beta^2\mu$ has been set to zero, using that $\tilde{\mu}$ is very close to $\mu$ in this temperature regime as discussed above (we also assume that $t$ is small, as justified below).

Following the analysis of the harmonic oscillator at zero temperature, we then introduce the scaled dimensionless vectors $\mathbf{a}$ and $\mathbf{b}$ defined via $\mathbf{a}' = w_N \mathbf{a}$ and $\mathbf{b}' = w_N \mathbf{b}$, where the width $w_N$ has the dimension of length, and is determined as follows. We impose that both the second [of order $O(t)$] and third term [of order $O(t^3)$] in Eq. (280) are of order unity. This determines both
$w_N$ and the typical time scale $t_N$ as

$$w_N = \frac{|\nabla V(\mathbf{r}_{\text{edge}})|^{-1/3} \hbar^{\frac{2}{3}}}{(2m)^{\frac{1}{3}}}, \quad t_N = |\nabla V(\mathbf{r}_{\text{edge}})|^{-2/3} (8m\hbar)^{\frac{1}{3}}, \qquad (281)$$

where the amplitudes have been chosen for later convenience. With this choice it is clear that the diffusion term $\frac{m}{2\hbar t}(\mathbf{a}' - \mathbf{b}')^2 \sim w_N^2/t_N \sim O(1)$ is also of order unity. Finally, for the term containing temperature to be also of $O(1)$ and for the parameter $b$ to be also of order unity, we must scale the temperature as a function of $N$ as follows:

$$\beta = \frac{bt_N}{2^{\frac{2}{3}}\hbar} = \frac{b(2m)^{\frac{1}{3}}}{\hbar^{\frac{2}{3}}|\nabla V(\mathbf{r}_{\text{edge}})|^{\frac{2}{3}}}. \qquad (282)$$

Specializing to the harmonic oscillator $V(r) = \frac{1}{2}m\omega^2 r^2$, we obtain

$$\beta = \frac{b(2m)^{\frac{1}{3}}}{\hbar^{\frac{2}{3}}(m\omega^2 r_{\text{edge}})^{2/3}} = \frac{b}{\hbar\omega} \frac{1}{[\Gamma(1+d)N]^{1/(3d)}}. \qquad (283)$$

In particular in $d=1$,

$$\beta = \frac{b}{\hbar\omega N^{1/3}} \qquad (284)$$

in agreement with the result of Sec. V B [see Eq. (122)]. In addition, the above formula (281) for the parameter $w_N$ when applied to the harmonic oscillator yields back the expression in Eq. (202) in any $d$. Note that for more general potentials $V(r) \sim r^p$ one finds from (282) that $\beta \sim b/T_{\text{edge}}$ where the temperature scale which controls the thermal fluctuations at the edge is

$$T_{\text{edge}} \sim N^{\frac{4(p-1)}{3d(p+2)}} \qquad (285)$$

which generalizes the scale $T_{\text{edge}} \sim N^{1/3}$ for $p=2$ and $d=1$.

Defining $\tau = t/t_N$ we can now rewrite

$$K_{\bar{\mu}}(\mathbf{x}, \mathbf{y}) \approx \frac{1}{C_d w_N^d} \int_{-\infty}^{\infty} \frac{du}{\exp(-bu) + 1} \int_{\Gamma} \frac{d\tau}{2\pi i} \frac{1}{\tau^{\frac{d}{2}}} \exp \left[ -\frac{(\mathbf{a} - \mathbf{b})^2}{2^{\frac{8}{3}}\tau} - \tau \left( \frac{a_n + b_n + 2u}{2^{\frac{1}{3}}} \right) + \frac{\tau^3}{3} \right], \qquad (286)$$

where $C_d = \pi^{d/2} 2^{(4d-2)/3}$. Using the integral representation of the diffusive propagator in Eq. (225) and the one of the Airy function in Eq. (206), we obtain the following scaling form for the edge kernel at finite temperature in arbitrary dimension $d$ for two points near the edge:

$$K_{\bar{\mu}}(\mathbf{x}, \mathbf{y}) \approx \frac{1}{w_N^d} K_{d,b}^{\text{edge}} \left( \frac{\mathbf{x} - \mathbf{r}_{\text{edge}}}{w_N}, \frac{\mathbf{y} - \mathbf{r}_{\text{edge}}}{w_N} \right), \qquad (287)$$

where we recall that $w_N = |\nabla V(\mathbf{r}_{\text{edge}})|^{-\frac{1}{3}} \hbar^{\frac{2}{3}} (2m)^{-\frac{1}{3}}$, and the scaling function is given by

$$K_{d,b}^{\text{edge}}(\mathbf{a}, \mathbf{b}) = 2^{\frac{2}{3}} \int_{-\infty}^{\infty} \frac{du}{\exp(-bu) + 1} \int \frac{d\mathbf{q}}{(2\pi)^d} e^{-i\mathbf{q}\cdot(\mathbf{a}-\mathbf{b})} \text{Ai}\left(2^{\frac{2}{3}}q^2 + \frac{a_n + b_n + 2u}{2^{1/3}}\right), \qquad (288)$$

which depends on a single dimensionless parameter $b = \hbar^{\frac{2}{3}} |\nabla V(\mathbf{r}_{\text{edge}})|^{\frac{2}{3}} (2m)^{-\frac{1}{3}}/T$ and the dimension of space $d$. Note that this result is independent of the precise shape of the potential (within the broad class to be discussed below) and is thus universal. The dependence on the potential enters only in the width parameter $w_N$, and the dimensionless inverse temperature $b$.

This representation can be expressed in an alternative form. We first split $\mathbf{q} = \mathbf{q}_t + q_n \mathbf{n}$ where $\mathbf{n} = \mathbf{r}_{\text{edge}}/r_{\text{edge}}$ is the direction normal to the edge, and we similarly split $\mathbf{a} = \mathbf{a}_t + a_n \mathbf{n}$ and $\mathbf{b} = \mathbf{b}_t + b_n \mathbf{n}$. We then carry out the integral over $q_n$ using the identity (231). We obtain the following alternative result for the scaling function of the edge kernel at finite temperature in space dimension $d$ as

$$K_{d,b}^{\text{edge}}(\mathbf{a}, \mathbf{b}) = \int_{-\infty}^{\infty} \frac{du}{\exp(-bu) + 1} \int \frac{d\mathbf{q}_t}{(2\pi)^{d-1}} e^{-i\mathbf{q}_t \cdot (\mathbf{a}_t - \mathbf{b}_t)} \text{Ai}(a_n + \mathbf{q}_t^2 + u) \text{Ai}(b_n + \mathbf{q}_t^2 + u) \qquad (289)$$

$$= \int \frac{d\mathbf{q}_t}{(2\pi)^{d-1}} e^{-i\mathbf{q}_t \cdot (\mathbf{a}_t - \mathbf{b}_t)} K_b^{\text{edge}}(a_n + \mathbf{q}_t^2, b_n + \mathbf{q}_t^2), \qquad (290)$$

where $K_b^{\text{edge}}$ is precisely the scaled edge kernel at finite temperature in $d=1$, given in Eq. (134).

*Average edge density at finite temperature.* The edge density is obtained simply by setting $\mathbf{a} = \mathbf{b}$ in Eq. (290). This leads to

$$N \rho_N(\mathbf{x}) = \frac{1}{w_N^d} F_{d,b}^{\text{edge}} \left( \frac{r - r_{\text{edge}}}{w_N} \right), \qquad (291)$$
where the scaling function in dimension $d$ can be obtained in terms of the one in $d=1$ as

$$F_{d,b}^{\text{edge}}(z) = \int \frac{d\mathbf{q}_i}{(2\pi)^{d-1}} F_{1,b}^{\text{edge}}(z + \mathbf{q}_i^2), \quad (292)$$

where the scaling function in $d=1$ is given in Eq. (126), which we recall here for convenience

$$F_{1,b}^{\text{edge}}(s) = \int_{-\infty}^{+\infty} du \frac{\text{Ai}(s+u)^2}{1+e^{-bu}}. \quad (293)$$

Note that for $d=1$ there is no integration over $\mathbf{q}_i$ to perform and (292) reduces to an identity. For $d > 1$ one can perform an angular integration to obtain

$$F_{d,b}^{\text{edge}}(z) = \frac{2^{2-d} \pi^{\frac{1-d}{2}}}{\Gamma(\frac{d-1}{2})} \int_{0}^{+\infty} dq \, q^{d-2} F_{1,b}^{\text{edge}}(z+q^2). \quad (294)$$

**Zero-temperature limit.** In the limit $T \to 0$, i.e., $b \to +\infty$, the Fermi factor converges to a Heaviside step function

$$\frac{1}{\exp(-bu) + 1} \to \theta(u) \quad (295)$$

and $\tilde{\mu} \to \mu$. The first representation of the finite-temperature kernel in Eq. (288) then recovers exactly the kernel obtained for the harmonic oscillator and given in Eq. (227). The second representation (289) leads to the alternative form for the zero-temperature edge kernel in terms of the Airy kernel (37) as

$$K_d^{\text{edge}}(\mathbf{a}, \mathbf{b}) = \int \frac{d\mathbf{q}_i}{(2\pi)^d} e^{-i\mathbf{q}_i \cdot (\mathbf{a}_i - \mathbf{b}_i)} K_{\text{Airy}}(a_n + \mathbf{q}_i^2, b_n + \mathbf{q}_i^2). \quad (296)$$

This thus generalizes to any $d$ the result given in Eq. (37) for $d=1$, i.e., the standard Airy kernel, and coincides, in an equivalent alternative form, with the harmonic oscillator result given in Eq. (227).

**Validity of the method.** We now return to the question on the range of validity of this universal edge kernel at finite temperature. In the derivation we essentially made two approximations: (i) a short-time expansion keeping terms only up to order $O(t^3)$, (ii) a gradient expansion of the potential, assuming that higher-order terms are subdominant for large $N$. To examine the validity of these two points, let us first focus on the potentials of the form $V(r) \sim r^p$ (with $p > 0$). In this case, we know from (253) that $r_{\text{edge}} \sim N^{\frac{2}{d(p+2)}}$. It follows from (281) since $|\nabla V| \sim r_{\text{edge}}^{p-1}$, that

$$w_N \sim N^{-\frac{2}{3d}(p-1)/(p+2)}, \quad t_N \sim w_N^2. \quad (297)$$

Furthermore, one can check that for $p > 0$ all the neglected terms are indeed subdominant for large $N$ (see Appendix A). For instance, the leading term $O(t^2)$, from Eq. (278), scales as $\sim \frac{t^2}{m} \nabla^2 V \sim t_N^2 r_{\text{edge}}^{p-2} \sim N^{-2/(3d)}$ which is negligible compared to the main $O(1)$ terms at large $N$. Similarly, it is easy to check that the neglected second term in $S_3$ in Eq. (245) is indeed small since $\frac{\hbar^2 \nabla^4 V}{m |\nabla V|^2} \sim r_{\text{edge}}^{p-2} \sim N^{-2/d}$. Finally, the gradient expansion is controlled by the parameter $w_N/r_{\text{edge}} \sim N^{-2/(3d)}$ which is small for all fixed $p$ at large $N$. Our conclusion is thus that all spherically symmetric polynomial potentials with

leading degree $p > 0$ do satisfy the validity criteria for our universal results at the edge.

We expect that the class of such potentials is actually much broader. The precise conditions are analyzed in Appendix A. However, there are well-identified cases where this universality breaks down, for instance, for wall-type potentials. One example is a box with an infinite hard wall. Another example is the potential $V(r) = \frac{1}{r^2} + r^2$ potential. At zero temperature, the latter is known to be related to Wishart matrices, which have different edge properties (close the origin) than the GUE ensemble. In this case, the limiting kernel is the so-called Bessel kernel [37,82].

In the case of the potential $V(r) \sim r^p$ with $0 < p < 1$, some additional peculiarities arise. First, one finds, from the above estimates, that the typical width of the edge region $w_N$ increases with increasing $N$. In addition, the temperature scale $T_{\text{edge}}$ defined in Eq. (285) decreases with increasing $N$ which is consistent with the fact that the potential is rather shallow so even a small temperature is sufficient to excite the system. Nevertheless, as discussed in Appendix A, the conditions for the universality class of the $p=2$ case studied in this paper seem to still hold, at least at a perturbative level. To assess more precisely the validity of this statement for $0 < p < 1$ would require further studies.

# VIII. DISCUSSION AND OPEN PROBLEMS

## A. Bosonization in the bulk and interactions: Beyond the sine kernel

In the bulk of the Fermi gas, the density $N\rho_N(x)$ varies very slowly compared to the typical spacing between fermions and other, more standard, methods can be applied. As discussed in the text, see around formula (195), the universal sine-kernel correlations (in $d=1$) and its high-dimensional generalizations can also be derived using the LDA method.

In $d=1$, another tool can be applied: the bosonization technique. We will briefly recall here the results which can be obtained from this method. One motivation is that it also allows to treat the case of interacting fermions, hence to address, to some extent, the question of universality in presence of interactions. Furthermore, although well known in the condensed matter community, these results do not seem widely known to mathematicians working on random matrices.

The bosonization method allows one to represent fermions in $d=1$ with uniform local density $\rho(x) := N\rho_N(x) \approx \rho_0$, as some “exponentials” of two conjugated bosonic fields, described by a quadratic Hamiltonian (see, e.g., [83] for a pedagogical introduction). This representation is exact for free fermions with a linear dispersion relation. For more general dispersion relations, and for interacting fermions, it remains an accurate effective description in the hydrodynamic limit, i.e., on spatial scales $(x-x')\rho_0 \gg 1$: this is the so-called Luttinger liquid (LL). The effective quadratic Hamiltonian is parametrized by the (renormalized) Luttinger parameter $K$, which contains all the information about the large scales. The special value $K=1$ corresponds to noninteracting fermions, while attraction leads to $K > 1$ and repulsion to $K < 1$. Using these methods it was shown [83,84] that the correlation
function of the density at $T = 0$ is given by [85]

$$ \langle \rho(x) \rho(0) \rangle_0 \approx \rho_0^2 \left[ 1 - \frac{2K}{(2\pi\rho_0 x)^2} + \sum_{m=1}^{+\infty} A_m (\rho_0 x)^{-2Km^2} \times \cos(2\pi m \rho_0 x) \right], \quad (298) $$

while the correlation function of the fermionic field (which is the analog of the kernel) is

$$ \langle \Psi^{\dagger}(x) \Psi(0) \rangle_0 \approx \rho_0 \sum_{m=0}^{+\infty} C_m (\rho_0 x)^{-\frac{1}{2K}-2K(m+\frac{1}{2})^2} \times \sin \left[ 2\pi \left( m + \frac{1}{2} \right) \rho_0 x \right]. \quad (299) $$

These formulas (298) and (299) are valid for $\rho_0 x \gtrsim 1$. Here, $A_1$ in Eq. (298) and $C_0$ in Eq. (299) represent the leading behaviors at large $\rho_0 x$, while the terms $A_m$, $m \ge 2$ and $C_m$, $m \ge 1$ represent the contributions of higher harmonics (often neglected in LL studies).

For noninteracting fermions, $K = 1$, $C_0 = 1$, and all $C_m = 0$ for $m \ge 1$, and the expression in Eq. (299) becomes exact. It is precisely the sine kernel $\sin(\pi\rho_0 x)/x = \ell(x)^{-1}K^{\text{bulk}}[x/\ell(x)]$ with $\ell(x) = 2/(\pi\rho_0)$ and $K^{\text{bulk}}(y) = \sin(2y)/\pi y$, proved in this paper to further hold for fermions in a trap (in the bulk), using the mapping to RMT. In presence of interactions we see that the correlation function of the fermionic field in Eq. (299) in the ground state now decays at large $x$ as

$$ \langle \Psi^{\dagger}(x) \Psi(0) \rangle_0 \sim \sin(\pi \rho_0 x) / x^{\eta}, \quad \eta = \frac{1}{2}(K + K^{-1}) \quad (300) $$

with a nonuniversal prefactor. The exponent $\eta$ is thus always larger than for noninteracting fermions, and its precise value depends on the Luttinger parameter $K$, hence on the strength of the interactions.

Unfortunately, the standard bosonization methods fail near the edge where $k_F = \pi\rho_0$ vanishes. We have found in this paper that for noninteracting fermions in a trap, the correlation function of the fermion field is described by the Airy kernel. Deriving that result using bosonization techniques seems at present out of reach. Nevertheless, it may be possible, at least qualitatively, to recover the leading asymptotics of the left tail of that kernel, i.e., the limit where both points enter into the bulk. Indeed, in the noninteracting case, we observe that by changing $\pi\rho_0 x \to \int_x^{x_{\text{edge}}} dy k_F(y)$ and inserting $k_F(y) \sim \sqrt{x_{\text{edge}} - y}$ in Eq. (299), we obtain the Airy function asymptotics at large negative arguments. Note that since the edge regime is diluted, it is possible that the effect of (at least short-range) interactions is less important than in the bulk, and the edge universality is more robust. This, however, remains to be studied in detail.

## B. Connection to the KPZ equation

In a recent paper [11] we have unveiled a remarkable connection between the problem of noninteracting fermions in a one-dimensional trap at finite temperature and the continuum $1d$ Kardar-Parisi-Zhang (KPZ) growth equation at finite time, with the so-called droplet initial condition (also called curved geometry). The KPZ equation [86] describes the stochastic time evolution of the height field $h(x,t)$ of an interface, at

point $x \in \mathbb{R}$ and time $t$,

$$ \partial_t h(x,t) = \nu \partial_x^2 h(x,t) + \frac{\lambda_0}{2} [\partial_x h(x,t)]^2 + \sqrt{D\eta} h(x,t), \quad (301) $$

where $\nu > 0$ is the coefficient of diffusive relaxation, $\lambda_0 > 0$ is the strength of the nonlinearity, and $\eta(x,t)$ is a centered Gaussian white noise with correlator $\langle \eta(x,t)\eta(x',t') \rangle = \delta(x-x')\delta(t-t')$. From now on, we express height and time in the natural units [87]

$$ t^* = 2(2\nu)^5/(D^2\lambda_0^4), \quad x^* = (2\nu)^3/(D\lambda_0^2), \quad h^* = 2\nu/\lambda_0. \quad (302) $$

Here, we start from the “narrow wedge” initial condition, $h(x,0) = -w|x| + \ln(w/2)$, with $w \gg 1$, which gives rise to a curved (or droplet) mean profile $\langle h(x,t) \rangle = \langle h(0,t) \rangle - \frac{x^2}{4t}$ as time evolves. Equivalently to (301) the Cole-Hopf transformed field $Z(x,t) = e^{h(x,t)}$ satisfies the stochastic heat equation (with multiplicative noise)

$$ \partial_t Z(x,t) = \partial_x^2 Z(x,t) + \sqrt{2\eta(x,t)} Z(x,t), \quad (303) $$

with initial condition $Z(x,t=0) = \delta(x)$. The continuum KPZ equation [88] is usually defined by the equation (303) with the Ito convention, implying that the first moment $\langle Z(x,t) \rangle = Z_0(x,t)$ where $Z_0(x,t) := \frac{1}{\sqrt{4\pi t}} e^{-x^2/(4t)}$ is the free diffusion propagator [90]. This is the definition we use here. In the natural units defined above, one can define the scaled height at position $x=0$:

$$ \tilde{h}(0,t) = \frac{h(0,t) + \frac{t}{12}}{t^{1/3}}. \quad (304) $$

The exact results of Refs. [29–32] can then be expressed as follows. Let us define the time-dependent generating function

$$ g_t(s) = \langle \exp(-t^{1/3}[\tilde{h}(0,t) - s]) \rangle \quad (305) $$

of the rescaled height at $x=0$. It is expressed for all time $t > 0$ as a Fredholm determinant:

$$ g_t(s) = \det (I - P_s K_t^{\text{KPZ}} P_s), \quad (306) $$

$$ K_t^{\text{KPZ}}(x,y) = \int_{-\infty}^{\infty} \frac{\text{Ai}(z+x)\text{Ai}(z+y)}{e^{-t^{1/3}z} + 1} dz, \quad (307) $$

where $P_s$ is the projector on the interval $[s, +\infty)$. We can now compare the finite-time kernel for KPZ, $K_t^{\text{KPZ}}$ in Eq. (307), with the finite-temperature kernel $K_b^{\text{edge}}$ for the fermions at the edge, given in Eq. (134), and see that they are identical provided one identifies

$$ b := \frac{\hbar\omega N^{1/3}}{T} \equiv t^{1/3}, \quad (308) $$

i.e., large time in KPZ corresponds to zero temperature of the fermions, and high temperature for the fermions to small time in KPZ [in original KPZ units the right-hand side should be replaced by $(t/t^*)^{1/3}$]. The correspondence appears more direct if one introduces an additional random variable $\gamma$, independent from $h(0,t)$ and distributed according to the Gumbel distribution $P(\gamma) = e^{-\gamma}e^{-\gamma'}$. Using that $\langle \theta(x - \gamma) \rangle_{\gamma} = e^{-e^{-\gamma}}$, we can rewrite the generating function (305) and obtain the
identity in law

$$
\lim_{N \to +\infty, \text{fixed } b} \xi := \frac{x_{\text{max}}(T) - x_{\text{edge}}}{w_N} \equiv_{\text{in law}} \frac{h(0,t) + \frac{t}{12} + \gamma}{t^{1/3}}, \quad (309)
$$

where the random variables on the left- and right-hand sides
have identical PDF. On the left-hand side, the limit of large
N is taken at fixed b, and the identity in law then holds for
arbitrary t and b related by (308). In the limit of large time
(respectively low T), the Gumbel variable can be ignored in the
right-hand side and one recovers that both the scaled variable
ξ = x_max(T) - x_edge / w_N and the scaled height h̃(0,t) tend to the Tracy-
Widom GUE distribution [of CDF F₂(s) given in Eq. (48)].

In the limit of small time one can use the small-t expansion
obtained in Ref. [30]. Using the formulas (11) and (12)
therein (and below the formula for the skewness) with the
correspondence ln z → h(0,t) − ln Z₀(0,t) + t/12 and λ³ =
t/4 = b³/4, one first obtains the first three cumulants for the
KPZ height field h(0,t) at small t. From this, using (309),
together with the cumulants of the Gumbel distribution
⟨γⁿ⟩ = (n − 1)!ζ(n) for n ≥ 2 and ⟨γ⟩ = γₐ, as well as
the independence of the two random variables, we obtain the
high-temperature expansion of the cumulants of the variable
ξ, i.e., the scaled position of the rightmost fermion. The
corresponding formulas are displayed in Sec. VC 2.

It is interesting to note that the relation (309) between
random variables can be “inverted” to express the PDF of
the KPZ field. Indeed, one can also write [29,30,32]

$$
\tilde{h}(0,t) = u - \frac{\gamma'}{t^{1/3}}, \quad (310)
$$

where $\gamma'$ is (yet another) Gumbel random variable, and $u$ is
a random variable, independent of $\gamma'$ and of time-dependent
PDF $p_t(u)$ obtained as

$$
p_t(u) = \mathrm{Det}[I - P_u(B_t - \mathrm{Ai}\mathrm{Ai}^\dagger)P_u] - \mathrm{Det}[I - P_u B_t P_u], \quad (311)
$$

$$
B_t(r, r') = P \int_{-\infty}^{+\infty} dv \frac{\mathrm{Ai}(r+v)\mathrm{Ai}(r'+v)}{1 - e^{-t^{1/3}v}} \quad (312)
$$

in terms of the kernel $B_t$ and of the projector $\text{Ai}\text{Ai}^\dagger(r,r') =$
$\text{Ai}(r)\text{Ai}(r')$ and where $P$ denotes the Cauchy principal value.
The large time expansion, i.e., the corrections around the TW
distribution (or low-temperature expansion for the fermions)
is studied in Appendix B by using equivalently (306) or (311).

Our results establish a precise connection between free
fermions at finite temperature and the KPZ equation (301)
at finite time *t*. At this stage, this connection exists only for the
droplet initial condition of KPZ and also for the one-point
distribution of *h*(0,*t*) [the one-point distribution of *h*(*x*,*t*)
being identical up to a shift by its average —*x*²/(4*t*).] Let
us mention two other examples where fermions appear in the
context of models in the KPZ universality class. One is the
polynuclear growth model, also related to a zero-temperature
model of directed paths [27,92]. This problem is related to
fermions on a 1*d* lattice in a linear potential, with a time-
dependent slope. A second example was studied recently in a
model of semidiscrete directed polymers [93]. It is at present
unclear whether the connection between fermions and the KPZ

equation hides a deeper correspondence, e.g., extending to
many point correlations or different initial conditions.

C. Conclusion

In this paper, we have developed a unified framework to
study the statistical mechanics of *N* noninteracting fermions
trapped by a confining potential, in any dimension *d* and at
any finite temperature *T*. The trapping potential gives rise to
an edge, i.e., a distance from the center of the trap above
which the average density vanishes. Consequently, when *N*
is large, there are two distinct scaling regimes in the local
correlations: the “bulk regime,” i.e., near the center of the
trap where the density of fermions varies smoothly, and the
“edge regime,” i.e., close to the edge where the density is
vanishing and where the (quantum and thermal) fluctuations
are thus large. In the “bulk regime,” our method recovers and
puts on firmer basis the results of the standard local density
(or Thomas-Fermi) approximation (LDA) [4,9]. On the other
hand, at the edge, where the LDA fails, we have obtained a
detailed description of the correlations. Indeed, we have shown
that, even at finite temperature *T* > 0, the system, in the limit
of a large number of fermions, is a determinantal point process
characterized by a kernel which depends on both the dimension
*d* and the temperature *T* [see Eq. (296)], which generalizes
the Airy kernel (37), which is a fundamental object in RMT.
Remarkably, we have shown, using a path-integral represen-
tation of this kernel, that is universal, i.e., independent of the
details of the trapping potential, for a wide class of spherically
symmetric potentials *V*(|*x*|) which behaves at large distance
as *V*(|*x*|) ~ |*x*|*p*, with *p* > 0. In the special case *d* = 1,
we have studied in detail the cumulative distribution function
of the position of the rightmost fermion *x*<sub>max</sub>(*T*) at *T* > 0
and we have shown that this CDF in Eq. (137) (i) generalizes
the well-known Tracy-Widom distribution, which describes
the fluctuations of *x*<sub>max</sub>(*T* = 0) and (ii) displays a remarkable
connection with the (1+1)-dimensional Kardar-Parisi-Zhang
equation in a droplet geometry (309). Furthermore, in the
case of a harmonic potential in any *d* we have shown that
all our results also hold in momentum space [in particular,
the Tracy-Widom distributions describes the fluctuations of
the largest momentum (56)]. It is an open question as to
how the latter property extends to nonharmonic traps.

Therefore, we believe that the results presented here open
the way for interesting bridges between the techniques of
many-body physics and of random matrix theory. We hope
that our results will stimulate further works at this interface.
In addition, one of the main outcomes of our paper is a precise
set of predictions for systems of noninteracting cold fermions
at finite temperature. It would be exciting if our theoretical
predictions could be verified in cold-atom experiments, for
example, using the state-of-the-art quantum microscopes
[6–8], or more conventional time-of-flight experiments.

ACKNOWLEDGMENTS

We would like to thank I. Bloch, T. Giamarchi, C. Salomon,
and G. Salomon for useful discussion. In addition, we would
like to thank the Kavli Institute for Theoretical Physics (KITP)
at Santa Barbara where part of this wok was done. This research
was supported in part by the National Science Foundation
under Grant No. NSF PHY11-25915.
APPENDIX A: GENERAL SHORT-TIME EXPANSION OF THE IMAGINARY-TIME SCHRÖDINGER PROPAGATOR

Here we derive, using probabilistic methods, the general short-time expansion for the imaginary-time propagator of a one-body Hamiltonian with an arbitrary soft potential. The semiclassical expansion of the path integral for the quantum mechanical propagator, an expansion in $\hbar$, has been extensively studied in the literature [94]. However the short-time expansion has received less attention [81]. As shown in the text, and detailed again in the next appendix, this expansion allows the calculation of the kernel of trapped noninteracting fermions in the bulk and at the edge. We first show the following formula for the propagator:

*Proposition.* The solution of the imaginary-time Schrödinger equation in Eq. (166) can be written, for an arbitrary potential $V(\mathbf{x})$, as

$$G(\mathbf{x}, \mathbf{y}; t) = \left(\frac{m}{2\pi\hbar t}\right)^{\frac{d}{2}} \exp\left[-\frac{m}{2\hbar}(\mathbf{x}-\mathbf{y})^2\right] \left\langle\exp\left\{-\frac{t}{\hbar}\int_0^1 du \, V[\mathbf{x}(1-u) + \mathbf{y}u + \sqrt{D_0t}\mathbf{B}_u]\right\}\right\rangle_{\mathbf{B}}, \quad (A1)$$

where $D_0 = \hbar/m$ is the diffusion coefficient and $\mathbf{B}_u = \{B_{iu}\}_{i=1}^d$ is a $d$-dimensional Brownian bridge on the interval $[0, 1]$, i.e., a Gaussian process with mean zero and correlation function

$$\langle B_{iu} B_{ju'} \rangle_{\mathbf{B}} = \delta_{ij} g(u, u'), \quad g(u, u') = \min(u, u') - uu', \quad (A2)$$

hence with $\mathbf{B}_0 = \mathbf{B}_1 = 0$.

*Proof.* To show this, we use the Feynman-Kac formula [95] to express the propagator as a path integral

$$G(\mathbf{x}, \mathbf{y}; t) = \int_{\mathbf{X}(0)=\mathbf{x}} d[\mathbf{X}] \delta(\mathbf{X}_t - \mathbf{y}) \exp \left\{ -\frac{1}{\hbar} \int_0^t ds \left[ \frac{1}{2} m \left( \frac{d\mathbf{X}_s}{ds} \right)^2 + V(\mathbf{X}_s) \right] \right\}, \quad (A3)$$

over all paths staring at $\mathbf{x}$ at time $t=0$, with a delta function weight so that only paths which end at $\mathbf{y}$ at time $t$ contribute. We now make a shift of variables in the path integral

$$\mathbf{X}_s = \mathbf{x} \frac{(t-s)}{t} + \mathbf{y} \frac{s}{t} + \mathbf{Z}_s, \quad (A4)$$

under which the expression for the propagator becomes

$$G(\mathbf{x}, \mathbf{y}; t) = \exp \left[ -\frac{m}{2\hbar} (\mathbf{x} - \mathbf{y})^2 \right] \left\langle \exp \left\{ -\frac{1}{\hbar} \int_0^t ds V \left[ \mathbf{x} \frac{(t-s)}{t} + \mathbf{y} \frac{s}{t} + \mathbf{Z}_t \right] \right\} \right\rangle_{\mathbf{Z}} \times \mathcal{N}_t, \quad (A5)$$

where $\mathcal{N}_t = \int_{\mathbf{Z}(0)=\mathbf{0}} d[\mathbf{Z}] \delta[\mathbf{Z}_t] \exp\left[-\frac{m}{2\hbar} \int_0^t ds \left(\frac{d\mathbf{Z}_s}{ds}\right)^2\right]$ is a normalization. In simplifying the kinetic energy contribution we have used that only paths where $\mathbf{Z}_t = \mathbf{0}$ contribute to the path integral, and we denote by $\langle \dots \rangle_{\mathbf{Z}}$ an average with respect to the normalized measure on the paths $\mathbf{Z}$ given by

$$P(\mathbf{Z}) = \frac{1}{\mathcal{N}_t} \delta[\mathbf{Z}_t] \exp \left[ -\frac{m}{2\hbar} \int_0^t ds \left( \frac{d\mathbf{Z}_s}{ds} \right)^2 \right]. \quad (A6)$$

To evaluate the normalization $\mathcal{N}_t$, we apply the Feynman-Kac formula again and see that the solution is given by the solution of the $d$-dimensional (free) diffusion equation, with a delta function initial condition. This gives $\mathcal{N}_t = (\frac{m}{2\pi\hbar t})^{\frac{d}{2}}$. The process $\mathbf{Z}$ with the measure given in Eq. (A6) is in fact a $d$-dimensional Brownian bridge, a Gaussian process with zero mean and temporal correlation function

$$\langle Z_{is} Z_{js'} \rangle_{\mathbf{Z}} = \delta_{is} D_0 \left[ \min(s, s') - \frac{ss'}{t} \right]. \quad (A7)$$

To show this, we compute the generating functional

$$g[\lambda] := \left\langle \exp \left[ -i \int_0^t \lambda(s) \cdot \mathbf{Z}_s \right] \right\rangle_{\mathbf{Z}} = \frac{\langle \delta(\mathbf{W}_t) \exp \left[ -i \int_0^t \lambda(s) \cdot \mathbf{W}_s \right] \rangle_{\mathbf{W}}}{\langle \delta(\mathbf{W}_t) \rangle_{\mathbf{W}}}, \quad (A8)$$

where $\langle \dots \rangle_{\mathbf{W}}$ denotes an average with respect to an unconstrained Brownian motion with correlation function

$$\langle W_{is} W_{js'} \rangle_{\mathbf{W}} = D_0 \min(s, s'). \quad (A9)$$

In the second equality in Eq. (A8) we enforced $\mathbf{W}_t = \mathbf{0}$ by inserting a delta function. The averaging over $\mathbf{W}$ can be carried out by using a Fourier representation of the delta function $\delta(\mathbf{W}_t) = \int \frac{d\mathbf{q}}{(2\pi)^d} \exp(-i\mathbf{q} \cdot \mathbf{W}_t)$. Inserting this representation, both in the denominator and in numerator of Eq. (A8), allows the Gaussian averaging to be carried out, yielding

$$g[\lambda] = \frac{\int \frac{d\mathbf{q}}{(2\pi)^d} \exp \left[ -\frac{D_0}{2} \int_0^t ds \, ds' \min(s, s') \lambda(s) \cdot \lambda(s') - \frac{D_0 t}{2} \mathbf{q}^2 - \int_0^t ds \, D_0 s \mathbf{q} \cdot \lambda(s) \right]}{\int \frac{d\mathbf{q}}{(2\pi)^d} \exp \left( -\frac{D_0 t}{2} \mathbf{q}^2 \right)}. \quad (A10)$$