# Statistics of the maximal distance and momentum in a trapped Fermi gas at low temperature

David S. Dean, Pierre Le Doussal, Satya. N. Majumdar, Gregory Schehr

► To cite this version:

David S. Dean, Pierre Le Doussal, Satya. N. Majumdar, Gregory Schehr. Statistics of the maximal distance and momentum in a trapped Fermi gas at low temperature. Journal of Statistical Mechanics: Theory and Experiment, 2017, 2017 (6), pp.063301 (1-39). 10.1088/1742-5468/aa6dda . hal-01558139

HAL Id: hal-01558139

https://hal.science/hal-01558139

Submitted on 7 Jul 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>
Statistics of the maximal distance and momentum in a trapped Fermi gas at low
temperature

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

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

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

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

We consider $N$ non-interacting fermions in an isotropic $d$-dimensional harmonic trap. We compute analytically the cumulative distribution of the maximal radial distance of the fermions from the trap center at zero temperature. While in $d = 1$ the limiting distribution (in the large $N$ limit), properly centered and scaled, converges to the squared Tracy-Widom distribution of the Gaussian Unitary Ensemble in Random Matrix Theory, we show that for all $d > 1$, the limiting distribution converges to the Gumbel law. These limiting forms turn out to be universal, i.e., independent of the details of the trapping potential for a large class of isotropic trapping potentials. We also study the position of the right-most fermion in a given direction in $d$ dimensions and, in the case of a harmonic trap, the maximum momentum, and show that they obey similar Gumbel statistics. Finally, we generalize these results to low but finite temperature.

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

Recent experimental advances in cold atom systems, in one or higher dimensions, have raised interesting new questions about these basic and fundamental systems [1, 2]. One such simple system corresponds to N non-interacting spin-less fermions in a d-dimensional harmonic trap. Near the center of the trap the system is just a free Fermi gas, which is a very well studied system in condensed matter physics and it can be analyzed using standard methods such as the Local Density Approximation (LDA) [3, 4]. However, the trap confines the Fermi gas in a limited region of space and thus introduces a sharp edge in the Fermi gas, in the limit of large N. Near this edge, the density is small and consequently quantum and thermal fluctuations are strong, leading to new emergent edge physics. The nature of this physics at the edges has been the subject of a number of recent studies [5–11].

The system of interest here consists of N non-interacting spin-less fermions in an isotropic d-dimensional harmonic potential $V(\mathbf{x}_1, \cdots, \mathbf{x}_N) = \frac{1}{2}m\omega^2 \sum_{i=1}^N \mathbf{x}_i^2$, where $\mathbf{x}_i$ denotes the position of the i-th fermion in d dimensions. In the following, for convenience, we will set $m = \omega = \hbar = 1$. The main focus of this paper will be on quantum fluctuations, i.e., setting the temperature $T = 0$, however some preliminary results on the role of finite temperature will also be given. At $T = 0$, the system will be in its many-body quantum ground state with the ground-state wave function denoted by $\Psi_0(\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N)$. The quantum fluctuations are then characterized by the joint probability distribution function (JPDF), $P(\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N) = |\Psi_0(\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N)|^2$, normalized to unity.

In $d = 1$, this many-body ground state wave function can be computed explicitly [6, 7, 12]

$$P(x_1, x_2, \cdots, x_N) = |\Psi_0(x_1, x_2, \cdots, x_N)|^2 = \frac{1}{Z_N} \prod_{1 \le i < j \le N} (x_i - x_j)^2 e^{-\sum_{i=1}^N x_i^2}, \quad (1)$$

where $Z_N$ is a normalization constant. Thus in $1d$ the positions $x_i$'s behave exactly as the eigenvalues of a random matrix belonging to the Gaussian Unitary Ensemble (GUE). In this case, the average density of fermions (normalized to unity) defined as $\rho_N(x) = (1/N) \sum_{i=1}^N (\delta(x - x_i))_0$ (where $\langle \cdots \rangle_0$ refers to a ground state average with respect to the JPDF in Eq. (1)) converges in the large $N$ limit to the celebrated Wigner semicircular form [12–14]

$$\rho_N(x) = \frac{1}{\sqrt{N}} f_W \left( \frac{x}{\sqrt{N}} \right) , \quad f_W(x) = \frac{1}{\pi} \sqrt{2 - x^2}. \quad (2)$$

Thus, on an average, the fermions are confined in a finite interval of the real line, $[-x_{\text{edge}}, +x_{\text{edge}}]$ where $x_{\text{edge}} = \sqrt{2N}$. Another striking prediction of this mapping [8] is that the quantum fluctuations of the position of the rightmost fermion at $T = 0$, i.e., the probability distribution of $x_{\text{max}} = \max\{x_1, x_2, \cdots, x_N\}$, converges in the large $N$ limit to the Tracy-Widom (TW) distribution characterizing the distribution of the largest eigenvalue of a GUE random matrix [15].

In $d > 1$, while the direct connection to Random Matrix Theory (RMT) no longer holds, it has been shown recently that several methods from RMT can still be exploited to study the correlations of the Fermi gas in $d > 1$ [9, 11]. Notably, in the large $N$ limit, there is a determinantal structure in the problem which allows one to express any ground-state correlation function as a determinant involving a kernel, for all $d \ge 1$. This kernel has been analyzed in detail in the limit of large $N$ [11]. For instance, the average density of the fermions, which is spherically symmetric, converges in the large $N$ limit to

$$\rho_N(\mathbf{x}) = \rho_N(r) \approx \frac{1}{N(2\sqrt{\pi})^d \Gamma(\frac{d}{2} + 1)} (r_{\text{edge}}^2 - r^2)^{d/2} \theta(r_{\text{edge}} - r), \quad (3)$$

with

$$r_{\text{edge}} = 2^{1/2} [\Gamma(d+1)]^{1/2d} N^{1/2d} = \sqrt{2\mu}, \quad (4)$$

where $\mu$ denotes the Fermi energy. The average density in Eq. (3) is a generalization of the Wigner semicircular law in $d \ge 1$. Thus, like in $d = 1$, the Fermi gas is confined in a finite sphere of radius $r_{\text{edge}}$ around the trap center. As mentioned earlier, in $d = 1$, the observable $x_{\text{max}}$ (the position of the rightmost fermion) obeys the Tracy-Widom distribution in the large $N$ limit. It is then natural to ask what is the analogue of this observable in higher dimensions and what is its limiting distribution in the large $N$ limit?

A natural observable in higher dimensions, that generalizes the position of the rightmost fermion for $d = 1$, is the radial distance of the farthest fermion

$$r_{\text{max}} = \max(r_1, r_2, \cdots, r_N) \quad \text{where} \quad r_i^2 = \mathbf{x}_i \cdot \mathbf{x}_i. \quad (5)$$

In the ground state, with JPDF $P(\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N)$, $r_{\text{max}}$ is clearly a random variable and in this paper our main focus will be to study the cumulative probability

$$P(w, N) = \text{Prob.}(r_{\text{max}} \le w, N), \quad (6)$$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Sketch of one configuration of N non-interacting fermions in a two-dimensional isotropic harmonic trap, centered at the origin, at zero temperature. The shaded area indicates the average density profile $\rho_N(r)$ (3), which vanishes beyond $r_{edge}$ in the large N limit. The fermion positions are labeled by $\mathbf{x}_i$, with $i = 1, \dots, N$ and their JPDF is given by $P(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N) = |\Psi_0(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N)|^2$, with $\Psi_0(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N)$ the many-body ground state wave function. Here we focus on the fluctuations of the radial distance $r_{max}$ of the farthest fermion (5). In the large N limit, we show here that the PDF of $r_{max}$, properly centered and scaled converges to a Gumbel distribution [see Eq. (7)].</figcaption></figure>

in all dimensions. Although we will obtain a formal solution for arbitrary N, the most interesting question concerns the large N limit. The answer to this question is straightforward in $d = 1$ where, using the exact mapping to RMT mentioned above, it is easy to show that for large N the limiting form of $P(w, N)$ converges to the squared TW distribution (see the discussion later). What is the analogous limiting form in $d > 1$? The main purpose of this paper is to address this question and present an exact result for this limiting distribution for all $d > 1$.

Let us summarize our main results. In the limit of large N and $d > 1$, we show that this cumulative distribution converges to the limiting scaling form

$$P(w, N) \underset{N \to \infty}{\longrightarrow} G\left(\frac{w - A_N}{B_N}\right), \quad \text{where } G(z) = e^{-e^{-z}}, \qquad (7)$$

where the scale factors $A_N$ and $B_N$ can be computed explicitly for large N. To leading order for large N, we show that

$$A_N \approx r_{\text{edge}} + \frac{1}{\sqrt{2\mu^{1/6}}} a_{\mu} \quad \text{and} \quad B_N \approx \frac{1}{2\sqrt{2} a_{\mu} \mu^{1/6}}, \qquad (8)$$

where $r_{edge} = \sqrt{2\mu}$ denotes the edge of the Fermi gas, $\mu \approx (\Gamma(1+d)N)^{1/d}$ is the Fermi energy and $a_\mu$ is given by

$$a_{\mu} \approx \left( \frac{d-1}{2} \ln \mu \right)^{2/3}. \qquad (9)$$

The limiting distribution $G(z)$ in Eq. (7) is a Gumbel distribution. In addition, we also study the distribution of $x_{max}$, the position of the rightmost fermion in $d$ dimensions, and again find a limiting Gumbel law. We also investigate the effect of finite temperature and show again the emergence of a Gumbel distribution for these observables ($r_{max}$ and $x_{max}$) although with different parameters.

It is well known that the Gumbel distribution emerges as a limiting distribution in the classical theory of extreme value statistics of independent and identically distributed (i.i.d.) random variables [16, 17]. In contrast, in the present case of fermions in a $d$-dimensional confining potential, the positions of the fermions are strongly correlated due to the Pauli exclusion principle. Thus it appears as a puzzle as to why the same Gumbel distribution for i.i.d. variables appears in this case. By decomposing the wave function into radial and angular coordinates, we will indeed unveil a “decorrelation” mechanism which explains the emergence of the Gumbel distribution, even in this strongly correlated system. It is interesting to note that a similar Gumbel limiting form also appears in the context of Ginibre matrices [12, 13, 18, 19], which are matrices with i.i.d. Gaussian entries (real, complex or quaternionic), but without any special symmetry of the matrix. In this case, the eigenvalues are generally distributed over the complex plane and
the distribution of the farthest eigenvalue from the origin, properly centered and scaled, converges for large N to a Gumbel distribution [20, 21] (see also [22] for a study of the large deviations beyond the Gumbel). For complex Ginibre matrices, the universality of the Gumbel law (with respect to different confining potentials), was established in Ref. [23]. The “decorrelating” mechanism for the Ginibre matrices thus appears to be similar to that for the fermion problem studied here.

One interesting property of the harmonic oscillator is that the position **x** and momentum **p** play completely symmetric roles. Hence, although we will not develop it here, all of our results for $r_{max}$ also carry out to $p_{max}$, the momentum of the fermion with the maximal momentum

$$p_{\max} = \max(p_1, p_2, \dots, p_N) \quad \text{where} \quad p_i^2 = \mathbf{p}_i \cdot \mathbf{p}_i. \qquad (10)$$

Hence we predict that this quantity will be governed by a Gumbel distribution in $d > 1$ for large N (while it involves the Tracy-Widom distribution for $d = 1$).

The paper is organized as follows. In section II, we recapitulate the general determinantal properties of non-interacting fermions in a harmonic potential and present the first derivation of the limiting distribution of the radial distance of the farthest fermion in $d > 1$, using the asymptotic properties of an underlying Fredholm determinant. In section III, we decompose the ground state wave function into radial and angular coordinates and present a second derivation of the same limiting distribution, that demonstrates the emergence of a “decorrelation” mechanism. In section IV, we discuss some preliminary extensions of our results and, in section V, we give our conclusions about this study.

## II. N FREE FERMIONS IN A d-DIMENSIONAL HARMONIC POTENTIAL: GROUND STATE

### A. Determinantal process

We consider N spin-less free fermions in a d-dimensional harmonic potential. The many-body Hamiltonian is given by

$$\hat{H} = \sum_{i=1}^{N} \hat{h}_i, \quad \text{where} \quad \hat{h}_i = -\frac{1}{2} \nabla_{\mathbf{x}_i}^2 + \frac{1}{2} r_i^2, \qquad (11)$$

where $\mathbf{x}_i$ is a d-dimensional vector denoting the position of the i-th fermion and $r_i^2 = \mathbf{x}_i \cdot \mathbf{x}_i$. For convenience, we have set $m = \hbar = 1$. The many-body eigenfunction $\Psi_E(\mathbf{x}_1, \dots, \mathbf{x}_N)$ of energy E satisfies the Schrödinger equation

$$\hat{H} \Psi_E(\mathbf{x}_1, \dots, \mathbf{x}_N) = E \Psi_E(\mathbf{x}_1, \dots, \mathbf{x}_N). \qquad (12)$$

Since the many-body Hamiltonian does not contain any interaction term, the many-body wavefunction can be expressed as a product of single particle eigenfunctions. The fermionic constraint however does not allow more than one particle in each single particle state. Consequently, the ground state wavefunction is given by the Slater determinant constructed from the N lowest single particle eigenfunctions

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

where $\mathbf{k}_i$ labels the single particle eigenfunction. In the ground-state, we fill up the single particle levels up to the Fermi energy which will be denoted by $\mu$. For the harmonic potential in d dimensions, the Fermi energy $\mu$ can be computed for large N and is given by [4, 9, 11]

$$\mu \approx [\Gamma(1+d)N]^{1/d}. \qquad (14)$$

Knowing exactly the many-body ground-state wavefunction $\Psi_0(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N)$ given in Eq. (13), the quantum joint probability density of the positions of the fermions in the ground state is given by

$$P(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N) = |\Psi_0(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N)|^2 = \frac{1}{N!} \det_{1 \le i,j \le N} [\psi_{\mathbf{k}_i}^*(\mathbf{x}_j)] \det_{1 \le i,j \le N} [\psi_{\mathbf{k}_j}(\mathbf{x}_i)]. \qquad (15)$$

Using $\det(A^T) \det(B) = \det(AB)$, the joint PDF can be re-expressed as a single determinant

$$P(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N) = \frac{1}{N!} \det_{1 \le i,j \le N} K_{\mu}(\mathbf{x}_i, \mathbf{x}_j), \qquad (16)$$
where we have defined the kernel $K_μ(\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}). \qquad (17)$$

This kernel plays a central role for the calculation of the correlations. For instance, the n-point correlation function $R_n(\mathbf{x}_1, \cdots, \mathbf{x}_n)$ defined as

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

can be expressed as an $n \times n$ determinant

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

This fact demonstrates that the free fermions in the ground state constitute a determinantal process [24, 25] in any dimension.

A simple consequence of this determinantal structure is that the average local density of fermions can also be expressed very simply in terms of the kernel $K_μ(\mathbf{x}, \mathbf{x})$ at identical points. Indeed, the average local density is defined as

$$\rho_N(\mathbf{x}) = \frac{1}{N} \left\langle \sum_{i=1}^{N} \delta(\mathbf{x} - \mathbf{x}_i) \right\rangle_0 \qquad (20)$$

where $\langle \cdots \rangle_0$ denotes the average w.r.t. the ground state quantum probability in Eq. (13). It is easy to show that

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

For isotropic quantum potentials, such as the harmonic potential considered here, the density $\rho_N(\mathbf{x}) = K_μ(\mathbf{x}, \mathbf{x})/N$, or equivalently $K_μ(\mathbf{x}, \mathbf{x})$ only depends on the radial coordinate $r = |\mathbf{x}|$. Hence

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

In the limit of large $N$, the average density is isotropic and has a simple form given by [9, 11]

$$\rho_N(r) \approx \frac{1}{N(2\sqrt{\pi})^d \Gamma(\frac{d}{2} + 1)} (r_{\text{edge}}^2 - r^2)^{d/2} \theta(r_{\text{edge}} - r), \qquad (23)$$

where $r_{\text{edge}}$ is given by [9, 11]

$$r_{\text{edge}} = 2^{1/2} \left[ \Gamma(d+1) \right]^{\frac{1}{2d}} N^{\frac{1}{2d}} = \sqrt{2\mu}. \qquad (24)$$

This bulk average density is well known and can also be derived by the standard semi-classical approximation (Thomas-Fermi or Local Density Approximation). However, this bulk form is only valid sufficiently far inside the edge at $r = r_{\text{edge}}$. It was pointed out recently that if one zooms in at the edge regime, the density profile does not vanish sharply as in Eq. (23), but instead gets smeared out over a scale $w_N$, due to finite $N$ corrections [9, 11]. The width of these edge fluctuations was computed explicitly in [9, 11]

$$w_N = \frac{1}{\sqrt{2}} [\Gamma(d+1) N]^{-\frac{1}{6d}} = \frac{1}{\sqrt{2}} \mu^{-\frac{1}{6}}. \qquad (25)$$

On this scale, i.e., when $r - r_{\text{edge}} \sim w_N$, the smeared density profile is given by $\rho_N(r) = K_μ(r, r)/N$ where [9, 11]

$$K_\mu(r, r) = \frac{1}{w_N^d} F_d \left( \frac{r - r_{\text{edge}}}{w_N} \right). \qquad (26)$$

The scaling function $F_d(s)$ was computed explicitly for all $d$ in Ref. [9, 11]. In this paper, we will need only the large $s$ asymptotics of $F_d(s)$ given by

$$F_d(s) \approx (8\pi)^{-\frac{d+1}{2}} s^{-\frac{d+3}{4}} \exp\left(-\frac{4}{3}s^{3/2}\right), \text{ as } s \to \infty. \qquad (27)$$
In this paper, we are interested in the statistics of the position of the farthest fermion from the trap center, in the ground state. More precisely, we define

$$r_{\max} = \max(r_1, r_2, \dots, r_N) \quad \text{where} \quad r_i^2 = \mathbf{x}_i \cdot \mathbf{x}_i, \qquad (28)$$

as the radial distance of the farthest fermion. This observable $r_{\max}$ is a random variable in the ground-state, due to the quantum fluctuations encoded in the joint PDF (15). Our goal is to compute the probability distribution function (PDF) of $r_{\max}$. It is however convenient first to consider the associated cumulative distribution function (CDF)

$$P(w, N) = \mathrm{Prob.}(r_{\max} \le w, N) = \mathrm{Prob.}[r_1 \le w, r_2 \le w, \dots, r_N \le w]. \qquad (29)$$

The PDF of $r_{\max}$ is thus given by $p(w, N) = \partial_w P(w, N)$. To compute $P(w, N)$, it is useful to introduce the indicator function $I_w(\mathbf{x}_i) = 1$ if $r_i \ge w$ and 0 otherwise, for each fermion. We can then express $P(w, N)$ as

$$P(w, N) = \left\langle \prod_{i=1}^{N} (1 - I_w(\mathbf{x}_i)) \right\rangle_0, \qquad (30)$$

where $\langle \dots \rangle_0$ refers to a ground state average with respect to the joint PDF in Eq. (15). One can evaluate this average on the right hand side (rhs) of Eq. (30) using the following identity valid for any determinantal process

$$\left\langle \prod_{i=1}^{N} (1 - f(\mathbf{x}_i)) \right\rangle_0 = \mathrm{Det}(I - f K_{\mu}), \qquad (31)$$

where $f(\mathbf{x})$ is an arbitrary function and Det denotes a Fredholm determinant, which is interpreted as

$$\mathrm{Det}(I - f K_{\mu}) = \exp[\ln \mathrm{Det}(I - f K_{\mu})] = \exp \left[ - \sum_{n=1}^{\infty} \frac{1}{n} \mathrm{Tr}[(f K_{\mu})^n] \right]. \qquad (32)$$

In this trace expansion, the term $\mathrm{Tr}[(f K)^n]$ is the $n$-dimensional integral

$$\mathrm{Tr}[(f K)^n] = \int d\mathbf{x}_1 \cdots \int d\mathbf{x}_n f(\mathbf{x}_1) K_\mu(\mathbf{x}_1, \mathbf{x}_2) f(\mathbf{x}_2) K_\mu(\mathbf{x}_2, \mathbf{x}_3) \cdots f(\mathbf{x}_n) K_\mu(\mathbf{x}_n, \mathbf{x}_1) \qquad (33)$$

where $K_\mu(\mathbf{x}, \mathbf{y})$ is the kernel defined in (17). Thus, choosing $f(\mathbf{x}) = I_w(\mathbf{x})$ in Eq. (31) gives the identity

$$P(w, N) = \mathrm{Det}(I - I_w K_{\mu}) = \exp \left[ - \sum_{n=1}^{\infty} \frac{1}{n} \mathrm{Tr}[(I_w K_{\mu})^n] \right] \qquad (34)$$

with

$$\mathrm{Tr}[(I_w K_{\mu})^n] = \int_{\mathcal{D}} d\mathbf{x}_1 \cdots d\mathbf{x}_n K_{\mu}(\mathbf{x}_1, \mathbf{x}_2) K_{\mu}(\mathbf{x}_2, \mathbf{x}_3) \cdots K_{\mu}(\mathbf{x}_n, \mathbf{x}_1) \qquad (35)$$

where the integration domain $\mathcal{D}$ corresponds to $r_i \ge w$, for each $i = 1, 2, \dots, N$. The reader should note that, up to this point, the formulas are valid for arbitrary finite $N$. In the next section we will use Eq. (34) to analyze the limiting form of $P(w, N)$ for large $N$.

## B. Limiting distribution of $r_{\max}$ for large $N$ in $d \ge 2$

Our goal is to compute the limiting form of $P(w, N)$ for large $N$ in all dimensions. As discussed in the introduction (see also later in section III), the case $d = 1$ is special, where the limiting distribution in the large $N$ limit, suitably centered and scaled, is given by a product of two Tracy-Widom distributions (see section III for details). In this subsection, we focus on $d \ge 2$ and we will see that the results here are very different from the $d = 1$ case.

Our starting point is the expression for $P(w, N)$ in Eq. (34). For large $N$, we anticipate that the PDF $p(w, N) = \partial_w P(w, N)$ must be peaked around the average value of $r_{\max}$. This average value, to leading order for large $N$, coincides with the edge of the density $r_{\mathrm{edge}}$ in Eqs. (23) and (24). Since $r_{\mathrm{edge}}$ is large for large $N$, we need to analyze the expression for $P(w, N)$ in Eq. (34) when $N$ and $w$ are both large. For large $w$ and $d \ge 2$, it turns out that $P(w, N)$ in Eq. (34) is dominated by the $n = 1$ term of the trace expansion – the proof of this fact is given in Appendix A. Keeping only this $n = 1$ term in (34) gives

$$P(w, N) \approx \exp[-\mathrm{Tr}(I_w K_{\mu})] = \exp\left[-\int_{\mathcal{D}} K_{\mu}(\mathbf{x}, \mathbf{x}) d\mathbf{x}\right] = \exp\left[-S_d \int_{w}^{\infty} K_{\mu}(r, r) r^{d-1} dr\right], \qquad (36)$$
where $S_d = 2\pi^{d/2}/\Gamma[d/2]$ is the surface area of a unit sphere in $d$ dimensions. Since $w$ is close to $r_{\text{edge}}$, we can use the scaling form in Eq. (26) for $K_{\mu}(r, r)$. Using Eq. (26) then gives

$$ P(w, N) \approx \exp \left[ -S_d \left( \frac{r_{\text{edge}}}{w_N} \right)^{d-1} \int_{\tilde{w}}^{\infty} F_d(s) ds \right], \qquad (37) $$

where

$$ \tilde{w} = \frac{w - r_{\text{edge}}}{w_N} = (w - \sqrt{2\mu}) \sqrt{2} \mu^{1/6}. \qquad (38) $$

Note that in Eq. (38) we have used the expressions for $r_{\text{edge}}$ and $w_N$, in terms of $\mu$, given respectively in Eq. (24) and (25). Using further that $r_{\text{edge}}/w_N = 2\mu^{2/3}$ in Eq. (152) gives

$$ P(w, N) \approx \exp \left[ -S_d 2^{d-1} \mu^{2(d-1)/3} \int_{\tilde{w}}^{\infty} F_d(s) ds \right]. \qquad (39) $$

For large $\mu$ and $d \ge 2$, the term $\mu^{2(d-1)/3}$ is large and consequently $P(w, N)$ is non zero only when the integral $\int_{\tilde{w}}^{\infty} F_d(s) ds$ is small and of order $\mu^{-2(d-1)/3}$. For this integral to be small, $\tilde{w}$ must be large. Hence we can evaluate this integral using the asymptotic form of $F_d(s)$ given in Eq. (27). Keeping only the leading term gives

$$ I(\tilde{w}) = \int_{\tilde{w}}^{\infty} F_d(s) ds \approx \frac{1}{2(8\pi)^{(d+1)/2}} \frac{1}{\tilde{w}^{(d+5)/4}} e^{-\frac{4}{3} \tilde{w}^{3/2}}. \qquad (40) $$

Substituting this result in Eq. (39) and exponentiating the prefactors inside the exponential, we get

$$ P(w, N) \approx \exp \left[ -B_d e^{-\frac{4}{3} \tilde{w}^{3/2} + \frac{2(d-1)}{3} \ln \mu - \frac{(d+5)}{4} \ln(\tilde{w})} \right], \qquad (41) $$

where $B_d$ is given by

$$ B_d = S_d \frac{2^{d-2}}{(8\pi)^{(d+1)/2}} = \frac{2^{-\frac{d+5}{2}}}{\Gamma(d/2)\sqrt{\pi}}. \qquad (42) $$

To obtain the limiting form of $P(w, N)$ in Eq. (41), we need to further center and suitably scale the variable $\tilde{w}$. For this purpose, we set

$$ \tilde{w} = a_{\mu} + b_{\mu} z \qquad (43) $$

where the two unknowns $a_{\mu}$ (the centering parameter or location) and $b_{\mu}$ (the scale parameter) will be chosen such that $P(w, N)$ becomes independent of $\mu$ (for large $\mu$) and is only a function of the variable $z$. We further anticipate (to be checked a posteriori) that, for large $\mu$, $a_{\mu}$ is large and $b_{\mu}$ is small, while $z$ is fixed ($\sim \mathcal{O}(1)$). Substituting this relation (43) in Eq. (41), expanding for large $\mu$, and keeping leading order terms, one gets

$$ P(w, N) \approx \exp \left[ -B_d e^{-\frac{4}{3} a_{\mu}^{3/2} + \frac{2(d-1)}{3} \ln \mu - \frac{(d+5)}{4} \ln(a_{\mu}) - 2 \sqrt{a_{\mu}} b_{\mu} z \right]. \qquad (44) $$

Hence we choose

$$ b_{\mu} = \frac{1}{2\sqrt{a_{\mu}}} \qquad (45) $$

and $a_{\mu}$ to be such that

$$ B_d e^{-\frac{4}{3} a_{\mu}^{3/2} + \frac{2(d-1)}{3} \ln \mu - \frac{(d+5)}{4} \ln(a_{\mu})} = 1. \qquad (46) $$

This then gives the limiting $N$-independent distribution

$$ P(w, N) \approx \exp[-e^{-z}]. \qquad (47) $$

Using (38), the result in Eq. (47) translates to

$$ \tilde{w} = (w - \sqrt{2\mu}) \sqrt{2} \mu^{1/6} = a_{\mu} + \frac{1}{2\sqrt{a_{\mu}}} z, \qquad (48) $$
where the random variable $z$ has the Gumbel CDF as in (47). Hence, in the large $N$ limit, one can express the original random variable $r_{\max}$ as an identity in law

$$r_{\max} = \sqrt{2\mu} + \frac{1}{\mu^{1/6}} \left[ \frac{1}{\sqrt{2}} a_{\mu} + \frac{1}{2\sqrt{2} a_{\mu}} Z_{G} \right], \quad (49)$$

where Prob.[$Z_G \le z$] = $\exp[-e^{-z}]$. The parameter $a_\mu$ is determined from Eq. (46) in the limit of large $\mu$. Solving Eq. (46) to leading orders for large $\mu$, we get

$$a_{\mu} \approx \left( \frac{d-1}{2} \ln \mu \right)^{\frac{2}{3}} + \frac{1}{2^{2/3} \left[ (d-1) \ln \mu \right]^{1/3}} \ln \left[ \frac{B_d}{\left( \frac{d-1}{2} \ln \mu \right)^{\frac{d+5}{6}}} \right], \quad (50)$$

where $B_d$ is given in Eq. (42). Equation (49), together with (50), is the main new result of this paper. Clearly, this Gumbel form for the limiting distribution holds strictly for $d > 1$.

### III. INTERPRETATION OF THE GUMBEL LAW

In the previous section we have investigated the statistics of $r_{\max} = \max\{r_1, r_2, \dots, r_N\}$ of $N$ non-interacting fermions in the ground state of a $d$-dimensional harmonic oscillator with $d > 1$. As discussed before, even though the fermions are non-interacting, the Pauli exclusion principle induces strong effective (repulsive) correlations between them. This makes the positions $r_i$'s of the fermions in the ground state strongly correlated. Therefore studying the statistics of $r_{\max}$ corresponds to studying the extreme value statistics of a strongly correlated set of random variables. For large $N$, we found that the limiting distribution of $r_{\max}$, suitably centered and scaled as in Eq. (49), is given by the Gumbel distribution. It is well known that the Gumbel law also appears in the context of classical extreme value statistics of *independent* and identically distributed (i.i.d.) random variables. So the fact that the Gumbel law also appears in a strongly correlated system of fermions in the ground state of a $d$-dimensional oscillator, with $d > 1$, is rather surprising. Is there any mechanism by which the effective random variables become “decorrelated” such that the Gumbel distribution emerges naturally? In the derivation presented in the previous section, this mechanism is not manifest.

In the present section, we present an alternative derivation of the same results as in Eqs. (49) and (50) using a completely different method that demonstrates explicitly how such a “decorrelating” mechanism emerges. The main idea is to decompose the ground state wave function in radial and angular sectors. For each angular quantum number, one effectively obtains a one-dimensional problem for a certain number of non-interacting fermions, in an effective quantum potential characterized by the angular quantum number, much like the hydrogen atom. In each of these angular sectors, the distribution of the farthest fermion position is non trivial. However, different angular sectors decouple and effectively one has to look at the maximum of a collection of independent, but non-identically distributed random variables. We will see below that indeed this mechanism eventually leads to the same Gumbel law, as in Eqs. (49) and (50).

To present this mechanism in a clear fashion, we first remind the reader how the single particle wavefunction can be decomposed in radial and angular sectors for a $d$-dimensional harmonic oscillator. Next we consider the $N$-particle problem and show that the CDF $P(w, N)$ of $r_{\max}$ factorizes into a product of determinants, each corresponding to a different angular quantum number. Each of these determinants can be mapped onto the distribution of the largest eigenvalue of a Laguerre-Wishart matrix. Finally, using the known results for the latter problem, we arrive at the result in Eqs. (49) and (50).

#### A. Single particle in a $d$-dimensional harmonic potential

We first recall the well known results for the quantum mechanics of a single particle in a $d$-dimensional harmonic potential. The Hamiltonian of the particle is given by:

$$\hat{h} = -\frac{1}{2} \nabla_{\mathbf{x}}^2 + \frac{1}{2} r^2, \quad (51)$$

where $\mathbf{x}$ is a $d$-dimensional vector denoting the position of the particle and $r^2 = \mathbf{x} \cdot \mathbf{x}$. The single particle eigenstates $\psi_{\mathbf{k}}(\mathbf{x})$ satisfy the Schrödinger equation

$$\hat{h} \psi_{\mathbf{k}}(\mathbf{x}) = \epsilon_{\mathbf{k}} \psi_{\mathbf{k}}(\mathbf{x}), \quad (52)$$
where $\epsilon_k$'s denote the single particle discrete energy levels. Note that we just use $k$ to label the single particle eigenstates (which are not plane waves). While it is easy to write down explicitly the eigenfunctions and eigenvalues in Cartesian coordinates, it is more convenient, for our purpose, to present the solution in spherical coordinates. The Hamiltonian operator can be split into a radial and an angular part. We denote the wave function in spherical coordinates by $\psi_k(r, \theta)$ where $\theta$ is a $(d-1)$-dimensional angular vector. In these coordinates the Schrödinger equation (52) becomes

$$ \frac{1}{2} \left[ -\frac{1}{r^{d-1}} \frac{\partial}{\partial r} \left( r^{d-1} \frac{\partial}{\partial r} \right) - \frac{\hat{\Delta}_{\theta}}{r^2} + r^2 \right] \psi_k(r, \theta) = \epsilon_k \psi_k(r, \theta), \qquad (53) $$

where $\hat{\Delta}_{\theta}$ is the Laplacian on the surface of a $d$-dimensional unit sphere [26]. Since the potential depends only on the radial coordinate $r = |x|$, we can write the solution using separation of variables as $\psi_k(r, \theta) = \phi(r)Y_L(\theta)$, where $L$ stands collectively for all the angular quantum numbers. For simplicity, we suppress, for the moment, the explicit quantum number dependence (i.e., $k$ dependence) of the eigenfunction $\phi(r)$. The angular part $Y_L(\theta)$ satisfies the eigenvalue equation

$$ \hat{\Delta}_{\theta} Y_{L}(\theta) = \lambda Y_{L}(\theta). \qquad (54) $$

It turns out that, while the eigenfunctions $Y_L(\theta)$ depend explicitly on all the angular quantum numbers denoted as $L$, the eigenvalues, in contrast, are labelled just by one single scalar quantum number [26], $\lambda = -l(l+d-2)$ where $l=0,1,2,...$. Each $l$-sector ($l>0$) is degenerate and the degeneracy is given by [26]

$$ g_d(l) = \frac{(2l+d-2)(l+d-3)!}{l!(d-2)!} = \frac{(2l+d-2)\Gamma(l+d-2)}{\Gamma(l+1)\Gamma(d-1)}. \qquad (55) $$

Note that $g_d(0) = 1$ in all dimensions $d \ge 1$. In $d=1$, it follows from Eq. (55) that

$$ g_1(0) = 1, \quad g_1(1) = 1, \quad g_1(l) = 0 \text{ for all } l \ge 2. \qquad (56) $$

This is consistent with the fact that in $d=1$, the quantum number $l$ has only two allowed values $l=0$ and $l=1$, corresponding to the angles $\theta=0$ and $\theta=\pi$ respectively.

These angular eigenfunctions in Eq. (54) satisfy the orthogonality condition

$$ \int d\theta Y_L(\theta) Y_{L'}(\theta) = \delta_{L,L'}. \qquad (57) $$

For a fixed $l$ the radial part satisfies the eigenvalue equation

$$ \left[ -\frac{1}{2} \frac{d^2}{dr^2} - \frac{d-1}{2r} \frac{d}{dr} + \frac{l(l+d-2)}{2r^2} + \frac{r^2}{2} \right] \phi(r) = \epsilon_k \phi(r). \qquad (58) $$

This equation for $\phi(r)$ (58) can be reduced to a standard one-dimensional Schrödinger equation by setting $\phi(r) = \chi(r)/r^{(d-1)/2}$ that eliminates the first order $d/dr$ term in Eq. (58) and we get

$$ -\frac{1}{2} \frac{d^2 \chi(r)}{dr^2} + \frac{1}{2} \left[ r^2 + \frac{(l + \frac{d-3}{2})(l + \frac{d-1}{2})}{r^2} \right] \chi(r) = \epsilon_k \chi(r) \quad , \quad r \ge 0. \qquad (59) $$

This $1d$ Schrödinger equation is solvable with discrete energy levels labelled by $n$ and $l$ [26] (see Fig. 2)

$$ \epsilon_k = \epsilon_{n,l} = 2n + l + \frac{d}{2}, \qquad (60) $$

where $n=0,1,2,...$. The associated radial eigenfunctions $\chi_{n,l}(r)$ are given explicitly by

$$ \chi_{n,l}(r) = A_{n,l} r^{l+\frac{d-1}{2}} L_n^{l+\frac{d-2}{2}}(r^2) e^{-\frac{r^2}{2}}, \qquad (61) $$

where $L_n^\alpha$ are generalized Laguerre polynomials (which are of degree $n$) and $A_{n,l}$ are normalization constants, such that the $\chi_{n,l}$ satisfy the orthogonality relation

$$ \int_0^\infty dr \, \chi_{n,l}(r) \chi_{n',l}(r) = \delta_{n,n'}. \qquad (62) $$
<figure><img src="image_3.png" /><figcaption>FIG. 2. Occupation of the energy levels in the plane (l, n) [see Eq. (63)] for N = 98 spin-less fermions in the ground state of the two-dimensional (d = 2) isotropic harmonic oscillator. The filled circles indicate occupied states while the open circles are unoccupied. In this case the energy levels are $\epsilon_{n,l} = 2n + l + 1$ and the degeneracy is $g_{n,l} = 2$ [see Eq. (55)], corresponding to $Y_L(\theta) = e^{\pm i l \theta}/\sqrt{2\pi}$, with $l = 0, 1, \dots$. Thus each dot actually corresponds to 2 distinct quantum states which implies that the Fermi energy (corresponding to the levels connected by the solid blue line) is $\mu = 13$. Note that the last energy level at the Fermi energy $\mu$ is fully occupied, the situation to which we restrict ourselves here (see text).</figcaption></figure>

Thus, finally, the eigenfunctions $\psi_k(\mathbf{x})$ of the d-dimensional harmonic oscillator, labelled by the quantum numbers $(n, \mathbf{L})$, are given by

$$ \psi_{n,l}(r, \theta) = \frac{1}{r^{\frac{d-1}{2}}} \chi_{n,l}(r) Y_{\mathbf{L}}(\theta), \qquad (63) $$

with an associated eigenvalue $\epsilon_{n,l} = (2n + l + d/2)$.

Note that in $d=1$, since only $l=0$ and $l=1$ are allowed, the additional potential term $l(l-1)/r^2$ in Eq. (59) vanishes and $\chi_{n,l}(r)$ satisfies the 1d Schrödinger equation in a harmonic potential, but restricted only to positive half-space ($r \ge 0$). The energy levels in Eq. (60) get divided into two sectors labeled by $l=0$ and $l=1$

$$ \epsilon_{n,0} = 2n + \frac{1}{2}, \quad \epsilon_{n,1} = (2n + 1) + \frac{1}{2}. \qquad (64) $$

One sees that the $l=0$ sector corresponds to the even eigenfunctions of the 1d harmonic oscillator in the full space, while the $l=1$ sector corresponds to the odd eigenfunctions.

## B. Many-body ground-state for N free fermions in a d-dimensional harmonic potential

*Ground-state wave function.* We now consider N spin-less free fermions in a d-dimensional harmonic potential. The many-body Hamiltonian is given by

$$ \hat{H} = \sum_{i=1}^{N} \hat{h}_{i}, \quad \text{where} \quad \hat{h}_{i} = -\frac{1}{2} \nabla_{\mathbf{x}_{i}}^{2} + \frac{1}{2} r_{i}^{2}, \qquad (65) $$

where $\mathbf{x}_i$ is a d-dimensional vector denoting the position of the i-th fermion and $r_i^2 = \mathbf{x}_i \cdot \mathbf{x}_i$. The many-body eigenfunction $\Psi_E(\mathbf{x}_1, \dots, \mathbf{x}_N)$ of energy $E$ satisfies the Schrödinger equation

$$ \hat{H} \Psi_E(\mathbf{x}_1, \dots, \mathbf{x}_N) = E \Psi_E(\mathbf{x}_1, \dots, \mathbf{x}_N). \qquad (66) $$

Each single particle Hamiltonian $\hat{h}_i$ has a discrete set of eigenvalues $\epsilon_{n_i, l_i} = (2n_i + l_i + d/2)$ as given above (60). Since the many-body Hamiltonian does not contain any interaction term, the many-body wavefunction can be expressed as a product of single particle eigenfunctions. The fermionic constraint however does not allow more than one particle in
each single particle state. Consequently, the ground state wavefunction is given by the Slater determinant constructed from the N lowest single particle eigenfunctions

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

where $\mathbf{k}_i = (n_i, \mathbf{L}_i)$ labels the single particle eigenfunction. In the ground-state, we fill up the single particle levels up to the Fermi energy which will be denoted by $\mu$. The ground-state energy $E_0$ is just the sum of the single particle energies of all filled levels.

Note that in $d > 1$ there is a large degeneracy of the single particle energy levels, i.e., many different couples $(n, l)$ have the same value of $\epsilon_{n,l}$, and in addition, for each such couple there is an angular degeneracy $g_d(l)$. In this work we will restrict ourselves to the case where the set of degenerate single particle energy levels with highest energy is fully occupied. In this case the N-body ground state is unique, and it is the case studied here (we refer the reader to Ref. [11] for further discussions of the degenerate ground states).

**Fermi energy.** Knowing explicitly the single particle energy levels $\epsilon_{n,l} = (2n+l+d/2)$, together with their degeneracies $g_d(l)$ given in Eq. (55), one can readily relate the number of fermions N to the Fermi energy $\mu$ as follows (see Fig. 2)

$$ N = \sum_{l, n \ge 0} g_d(l) \theta(\mu - \epsilon_{n,l}). \qquad (68) $$

Carrying out the sum over $n$, one obtains

$$ N = \sum_{l \ge 0} g_d(l) \left( 1 + \text{int} \left( \frac{\mu - l - d/2}{2} \right) \right) \theta \left( \mu - l - d/2 \ge 0 \right), \qquad (69) $$

where $\text{int}(x)$ denotes the integer part of $x$. We can easily estimate $\mu$ for large N. For this, one needs to distinguish the $d=1$ case and $d \ge 2$.

**Case $d=1$.** In $d=1$, the sum over $l$ in Eq. (69) runs over two terms, $l=0$ and $l=1$, with $g_1(0) = g_1(1) = 1$. It follows immediately from Eq. (69) that, for large N,

$$ \mu \approx N. \qquad (70) $$

**Case $d \ge 2$.** In the large N limit, $\mu$ is expected to be large. In this case, the sum on the right hand side of Eq. (69) is dominated by large $l$. For large $l$ and $d \ge 2$, using Stirling's formula in Eq. (55), the degeneracy $g_d(l)$ behaves asymptotically as,

$$ g_d(l) \approx \frac{2 l^{d-2}}{\Gamma(d-1)}. \qquad (71) $$

Substituting this behavior in the sum and rescaling $l = \mu\ell$, one can replace, for large $\mu$, the discrete sum by an integral over $\ell \in [0, 1]$. In this large $\mu$ limit, one can also replace $\text{int}((\mu - l - d/2)/2)$ by $\mu(1 - \ell)/2$. This gives

$$ N \approx \frac{\mu^d}{\Gamma(d-1)} \int_0^1 (1-\ell)^{d-2} dz = \frac{\mu^d}{\Gamma(d+1)}. \qquad (72) $$

One obtains

$$ \mu \approx (\Gamma(d+1) N)^{1/d}, \qquad (73) $$

thus recovering the result mentioned in (14). Interestingly, although the derivation holds only for $d > 1$, the result in Eq. (73) remains valid even for $d=1$, where it gives the correct result $\mu \approx N$ in (70).

Finally, we introduce another quantity that will be useful later. For a given $\mu$, it is clear from Eq. (68), that for fixed $l$ the sum over $n$ is cut off at the top due to the theta function (see also Fig. 2). Let $N_l(\mu)$ be this maximal allowed value. It is given by

$$ N_l(\mu) = \text{int} \left( \frac{\mu - l - d/2}{2} \right). \qquad (74) $$

For $d=1$, there are only two allowed values $l=0$ and $l=1$ and hence for large $\mu$, it follows from Eq. (74) that these two sectors have the same number of states

$$ N_{l=0}(\mu) \approx \frac{\mu}{2} \approx \frac{N}{2}, \quad N_{l=1}(\mu) \approx \frac{\mu}{2} \approx \frac{N}{2}. \qquad (75) $$
### C. The distribution of the position of the farthest fermion

We now turn to the CDF $P(w, N)$ of $r_{max}$ defined in Eq. (29). Using the representation in Eq. (30) and the joint PDF in Eq. (15), we get

$$P(w, N) = \left\langle \prod_{i=1}^{N} (1 - I_w(\mathbf{x}_i)) \right\rangle_0 \qquad (76)$$

$$= \frac{1}{N!} \int d\mathbf{x}_1 \dots \int d\mathbf{x}_N \det_{1\le i,j\le N} [\psi_{\mathbf{k}_i}^*(\mathbf{x}_j)] \det_{1\le i,j\le N} [\psi_{\mathbf{k}_i}(\mathbf{x}_j)] \prod_{i=1}^{N} (1 - I_w(\mathbf{x}_i)) . \qquad (77)$$

In section II, we had seen that $P(w, N)$ can be expressed as a Fredholm determinant [see Eq. (34)]. Here, we show that it can be expressed as a determinant of a finite matrix (as opposed to the functional or Fredholm determinant), using the powerful Andreief-Cauchy-Binet identity. According to this identity, one has

$$\int d\mathbf{x}_1 \dots \int d\mathbf{x}_N \prod_{i=1}^{N} h(\mathbf{x}_i) \det_{1\le i,j\le N} [f_i(\mathbf{x}_j)] \det_{1\le i,j\le N} [g_i(\mathbf{x}_j)] = N! \det_{1\le i,j\le N} \left[ \int d\mathbf{x} h(\mathbf{x}) f_i(\mathbf{x}) g_j(\mathbf{x}) \right], \qquad (78)$$

where $f_i(\mathbf{x})$, $g_j(\mathbf{x})$ and $h(\mathbf{x})$ are arbitrary integrable functions. Choosing $h(\mathbf{x}) = 1 - I_w(\mathbf{x})$ and setting $f_i = \psi_{\mathbf{k}_i}$ and $g_j = \psi_{\mathbf{k}_j}^*$, we immediately get

$$P(w, N) = \det_{1\le i,j\le N} \left[ \int d\mathbf{x} [1 - I_w(\mathbf{x})] \psi_{\mathbf{k}_i}^*(\mathbf{x}) \psi_{\mathbf{k}_j}(\mathbf{x}) \right] = \det_{1\le i,j\le N} \left[ \delta_{ij} - \int d\mathbf{x} I_w(\mathbf{x}) \psi_{\mathbf{k}_i}^*(\mathbf{x}) \psi_{\mathbf{k}_j}(\mathbf{x}) \right], \qquad (79)$$

where we have used the orthogonality of $\psi_{\mathbf{k}}$'s in the second equality. Using the parametrization of the quantum numbers $\mathbf{k} = (n, \mathbf{L})$, the expression above can be rewritten as

$$P(w, N) = \det_{1\le i,j\le N} \left[ \delta_{ij} - \int_{r\ge w} d\mathbf{x} \psi_{n_i, \mathbf{L}_i}^*(\mathbf{x}) \psi_{n_j, \mathbf{L}_j}(\mathbf{x}) \right], \quad \text{where } r^2 = \mathbf{x} \cdot \mathbf{x}, \qquad (80)$$

where $\psi_{n,\mathbf{L}}(\mathbf{x})$ is given in Eq. (63). Using $d\mathbf{x} = r^{d-1} dr d\theta$ in the integral on the rhs of Eq. (80) and the explicit form of $\psi_{n,\mathbf{L}}(\mathbf{x})$ in Eq. (63), we see that the factor $r^{d-1}$ cancels. Carrying out the angular integral and using the orthogonality condition in Eq. (57), one obtains

$$\int_{r \ge w} d\mathbf{x} \psi_{n_i, \mathbf{L}_i}^*(\mathbf{x}) \psi_{n_j, \mathbf{L}_j}(\mathbf{x}) = \delta_{\mathbf{L}_i, \mathbf{L}_j} \int_w^\infty dr \chi_{n_i, l_i}(r) \chi_{n_j, l_j}(r). \qquad (81)$$

Using this result in Eq. (81), we see that the determinant on the right hand side of Eq. (80) decouples into a product of determinants, each corresponding to a particular angular sector. Taking into account the degeneracy $g_d(l)$ in Eq. (55) associated to each value of the quantum number $l$, we get

$$P(w, N) = \prod_{l \ge 0}^{\text{int}(\mu - \frac{d}{2})} \left[ P_l(w, m_l) \right]^{g_d(l)}, \qquad (82)$$

where, for $d > 1$,

$$m_l = N_l(\mu) + 1 \quad \text{with} \quad N_l(\mu) = \text{int} \left( \frac{\mu - l - d/2}{2} \right), \qquad (83)$$

denotes the number of fermions in the states of angular momentum $l$. The product over $l$ in Eq. (82) runs up to $\mu - d/2$. By contrast, for $d = 1$, the product in Eq. (82) runs only over two values $l = 0$ and $l = 1$, with $g_1(0) = g_1(1) = 1$ and the corresponding values of $m_l$ are $N/2$ for each sector [see Eq. (75)].

In Eq. (82), $P_l(w, m_l)$ is an $m_l \times m_l$ determinant given by

$$P_l(w, m_l) = \det_{1\le i,j\le m_l} \left[ \delta_{i,j} - \int_w^\infty \chi_{n_i,l}(r) \chi_{n_j,l}(r) dr \right], \qquad (84)$$
where $n_i = 0, 1, \dots, m_l - 1$. In Eq. (84), the radial wavefunction $\chi_{n,l}(r)$ is given explicitly in Eq. (61). We will see below that this determinant in Eq. (84) can actually be interpreted as the cumulative distribution of the position of the rightmost fermion among a set of $m_l$ spin-less fermions in a one-dimensional quantum problem.

The effective 1d problem. In order to interpret $P_l(w, m_l)$ in (84), we consider an effective quantum 1d model on the positive real axis with a single particle Hamiltonian given by

$$ \hat{h}_{\text{eff}} = -\frac{1}{2} \frac{d^2}{dr^2} + \frac{1}{2} \left[ r^2 + \frac{\alpha(\alpha - 1)}{r^2} \right], \qquad (85) $$

with $\alpha = l + (d - 1)/2$. The single particle eigenfunctions are given by $\chi_{n,l}(r)$ given in Eq. (61) with eigenvalues $\epsilon_{n,l} = (2n + l + d/2)$ as in Eq. (60). We now consider $m_l = N_l(\mu) + 1$ non-interacting fermions with many-body Hamiltonian which is just the sum of the single particle Hamiltonians in Eq. (85). Indeed, this non-interacting fermion problem has already been studied in the context of non-intersecting fluctuating interfaces in 1+1 dimensions [27] and a relation to Laguerre-Wishart random matrices was found, as we recall below.

The many-body ground state wavefunction of this effective 1d fermion model can be written again as a Slater determinant of the first $m_l$ single particle eigenstates

$$ \Psi_{\alpha}(r_1, r_2, \dots, r_{m_l}) = \frac{1}{\sqrt{Z_{m_l, \alpha}}} \left[ \prod_{i=1}^{m_l} r_i^{\alpha} e^{-\frac{1}{2} r_i^2} \right] \left[ \det \prod_{1 \le j, k \le m_l} L_k^{\alpha-1/2}(r_j^2) \right], \qquad (86) $$

where $Z_{m_l, \alpha}$ is the normalization constant. The determinant of Laguerre polynomials can be expressed in terms of a Vandermonde determinant giving

$$ \left[ \det \prod_{1 \le j, k \le m_l} L_k^{\alpha-1/2}(r_j^2) \right] \propto \prod_{1 \le j < k \le m_l} (r_j^2 - r_k^2). \qquad (87) $$

The joint PDF of the positions of the $m_l$ Fermions in the ground state is given by

$$ p_{\alpha}(r_1, r_2, \dots, r_{m_l}) = |\Psi_{\alpha}(r_1, r_2, \dots, r_{m_l})|^2 = A_{m_l, \alpha} \left[ \prod_{i=1}^{m_l} r_i^{2\alpha} e^{-r_i^2} \right] \prod_{1 \le j, k \le m_l} (r_k^2 - r_j^2)^2, \qquad (88) $$

where $A_{m_l, \alpha}$ is a normalization constant, such that $\int dr_1 \dots \int dr_{m_l} p_{\alpha}(r_1, r_2, \dots, r_{m_l}) = 1$.

In this model, the CDF of the position of the rightmost fermion can be written as

$$ P_{\alpha, m_l}(w) = \text{Prob.}[\max\{r_1, r_2, \dots, r_{m_l}\} \le w] = \left\langle \prod_{i=1}^{m_l} (1 - I_w(r_i)) \right\rangle \qquad (89) $$

where $I_w(r_i) = 1$ if $r_i \ge w$ and 0 otherwise. The average $\langle \dots \rangle$ is carried out with respect to the joint distribution $p_{\alpha}(r_1, \dots, r_{m_l})$ in Eq. (88). Following similar steps as used from Eqs. (76) to (80), but specializing now to the 1d case, we get

$$ P_{\alpha, m_l}(w) = \det \left[ \delta_{i,j} - \int_w^\infty \chi_{n_i, l}(r) \chi_{n_j, l}(r) dr \right] = P_l(w, m_l), \qquad (90) $$

with $\alpha = l + (d - 1)/2$ and $m_l = N_l(\mu) + 1$. Thus the quantity $P_l(w, m_l)$ that appears in Eq. (82) can be interpreted as the CDF of the position of the rightmost fermion in the ground state of this effective 1d quantum model.

Making a change of variable $\lambda_i = r_i^2$, the joint PDF of the $\lambda_i$'s variables can be written from $p_{\alpha}(r_1, r_2, \dots, r_{m_l})$ in Eq. (88) as

$$ \rho_{\alpha}(\lambda_1, \lambda_2, \dots, \lambda_{m_l}) = \frac{1}{Z'_{m_l, l}} \left[ \prod_{i=1}^{m_l} \lambda_i^{\alpha-1/2} e^{-\lambda_i} \right] \prod_{1 \le j < k \le m_l} (\lambda_k - \lambda_j)^2 \qquad (91) $$

where $Z'_{m_l, l}$ is another normalization constant and we recall that $\alpha = l + (d - 1)/2$ and $m_l = N_l(\mu) + 1$. This is exactly the joint PDF of the eigenvalues of the Laguerre-Wishart random matrices with Dyson index $\beta = 2$ [12, 13].

Let us quickly recall the Laguerre-Wishart ensemble of random matrices. A Wishart matrix is an $m \times m$ product matrix of the form $W = X^\dagger X$ where $X$ is an $M \times m$ rectangular random matrix and $X^\dagger$ is its Hermitian conjugate,
which is of size $m \times M$ [28]. If one takes the distribution of the entries of $X$ to be Gaussian, drawn from the joint PDF

$$P(X) \propto \exp \left[ -\frac{\beta}{2} \mathrm{Tr}(X^{\dagger} X) \right], \qquad (92)$$

where the Dyson index $\beta = 1$ for real matrices and $\beta = 2$ for complex matrices. Without any loss of generality, we consider here $M \ge m$ (the case $M \le m$ can be analyzed by exchanging $M$ and $m$ [29]). The product matrix $W$ of size $m \times m$ has $m$ non-negative real eigenvalues $\lambda_1, \dots, \lambda_m$ whose joint PDF is given by

$$P_{\mathrm{LW}}(\lambda_1, \dots, \lambda_m) \propto \exp \left( -\frac{\beta}{2} \sum_{i=1}^{m} \lambda_i \right) \left[ \prod_{i=1}^{m} \lambda_i^{\frac{\beta}{2}(M-m+1)-1} \right] \prod_{1 \le j,k \le m} |\lambda_k - \lambda_j|^{\beta}, \qquad (93)$$

where the subscript LW refers to the Laguerre-Wishart ensemble. Even though in the original Wishart matrices $M$ and $m$ are integers, the more general Laguerre-Wishart ensemble corresponds to Eq. (93) where the parameter $M - m > -1$ can be any real number. We thus see that the distribution in Eq. (91) is the same as that of the eigenvalues of the $\beta = 2$ Laguerre-Wishart ensemble upon the identifications:

$$m = m_l = N_l(\mu) + 1 \quad \text{and} \quad M = M_l = m_l + l + \frac{d}{2} - 1, \qquad (94)$$

where $N_l(\mu)$ is given in Eq. (74). Going back to Eq. (90) and recalling that $\lambda_i = r_i^2$, we see that $P_l(w, m)$ can then be interpreted as the CDF of the maximal eigenvalue in the Laguerre-Wishart ensemble

$$P_l(w, m_l) = \mathrm{Prob.} \left[ \max \{ \lambda_1, \lambda_2, \dots, \lambda_{m_l} \} \le w^2 \right]. \qquad (95)$$

For later convenience, we define a random variable in each $l$-sector

$$r_{\max,l} = \sqrt{\max \{ \lambda_1, \lambda_2, \dots, \lambda_{m_l} \}}. \qquad (96)$$

Clearly, it follows from Eqs. (95) and (96) that

$$P_l(w, m_l) = \mathrm{Prob.}(r_{\max,l} \le w). \qquad (97)$$

Therefore from Eq. (82), we see that the CDF of $r_{\max}$, $P(w, N)$, is the product of CDFs of $(\mu - d/2)$ (interpreted as the integer part) independent random variables $r_{\max,l}$, one for each $l$. Each of these random variables is the square root of the maximum eigenvalue of a Laguerre-Wishart matrix with size parameters $m_l = N_l(\mu) + 1$ and $M_l = m_l + l + \frac{d}{2} - 1$. The Laguerre-Wishart ensembles labelled by $l$ are independent of each other. Thus $r_{\max}$ corresponds to the extreme of a set of independent but non-identically distributed random variables.

It is important to note that the above statements and formula are correct for arbitrary finite $N$ provided the last energy level is fully occupied, so that the ground state is unique. In the next subsection we study the asymptotics of $P(w, N)$ for large $N$, by analyzing the product structure in Eq. (82) and we will show that the same Gumbel law as in Eqs. (49) and (50) will emerge out of this analysis.

**Remark for $d = 1$:** We note that in $d = 1$, using the one to one correspondence between the fermion positions and the eigenvalues of the $N \times N$ GUE matrices, it follows that

$$P(w, N) = \mathrm{Prob.}[\max\{|x_1|, \dots, |x_N|\} \le w], \qquad (98)$$

where $x_i$'s denote the fermion positions (or equivalently the GUE eigenvalues). Going back to Eq. (82) and noting that, in $d = 1$, there are only two allowed values of $l$, namely $l = 0$ and $l = 1$ with degeneracies $g_1(0) = g_1(1) = 1$, it follows that the product in Eq. (82) has only two terms

$$P(w, N) = P_0(w, m_0) P_1(w, m_1) \qquad (99)$$

with $\mu = N - 1/2$, $m_0 = 1 + \mathrm{int}((N - 1)/2)$, $m_1 = 1 + \mathrm{int}((N - 2)/2)$. Thus for any $N$ there is a factorization into two sectors, which corresponds to two independent squared Laguerre ensembles [clarifying Eq. (96)]. This is reminiscent of the factorization of the joint PDF of the absolute values of the eigenvalues of GUE into two squared Laguerre ensembles (sometimes called chiral Laguerre ensembles) [32–34]. Again, these statements are valid for any finite $N$.
### D. Large N analysis of $P(w, N)$ in Eq. (82)

Our starting point is the expression for $P(w, N)$ in terms of a product, given in Eq. (82), where $P_l(w, m_l)$ is given in Eq. (95), which corresponds to the CDF of a Laguerre-Wishart matrices with size parameters $m_l = N_l(\mu) + 1$ and $M_l = m_l + l + \frac{d}{2} - 1$. In the following, in order to lighten the notations, we will drop the subscript $l$ from both $m_l$ and $M_l$, and denote them simply by $m$ and $M$.

Since the Fermi energy $\mu \propto N^{1/d}$ [see Eq. (73)], the limit of large $N$ corresponds to large $\mu$. Consequently, from Eq. (83), $m_l \equiv m$ is large for large $N$ and hence we need to analyze $P_l(w, m)$ for large $m$. The distribution of the largest eigenvalue $\lambda_{max}$ of Laguerre-Wishart matrices in the large $m$ limit has been well studied in the literature [30, 31]. Let us recall the main results that we will need here for our analysis of $P_l(w, m)$.

Consider a complex $(\beta = 2)$ Laguerre-Wishart ensemble with size parameters $m, M$ such that both $m$ and $M$ are large with their ratio held fixed. In this limit, the largest eigenvalue $\lambda_{max}$ approaches to an identity in law

$$ \lambda_{max} = A_{m,M} + \sigma_{m,M} \chi_2, \qquad (100) $$

where

$$ A_{m,M} = \left(\sqrt{m} + \sqrt{M}\right)^2 \quad \text{and} \quad \sigma_{m,M} = \left(\sqrt{m} + \sqrt{M}\right) \left(\frac{1}{\sqrt{m}} + \frac{1}{\sqrt{M}}\right)^{1/3} \qquad (101) $$

and $\chi_2$ is a random variable, independent of $m$ and $M$, whose CDF is given by the Tracy-Widom distribution

$$ \text{Prob.} (\chi_2 \le s) = \mathcal{F}_2(s). \qquad (102) $$

This function $\mathcal{F}_2(s)$ can be expressed in terms of a special solution of a Painlevé II equation and has the large $s \to \infty$ asymptotic tail

$$ \mathcal{F}_2(s) \sim 1 - \frac{1}{16\pi s^{\frac{3}{2}}} \exp\left(-\frac{4}{3}s^{\frac{3}{2}}\right) \quad \text{as} \quad s \to \infty. \qquad (103) $$

To adapt these results to our problem, we proceed as follows. Let us first express $A_{m,M}$ and $\sigma_{m,M}$ in terms of $\mu$ and $l$ for large $\mu$. From Eq. (83) for large $\mu$, we can ignore the integer part and write

$$ m \approx \frac{1}{2} \left( \mu - l - \frac{d}{2} + 2 \right) \approx \frac{1}{2} (\mu - l), \qquad (104) $$

anticipating that $\mu - l$ will be large (to be verified a posteriori). Similarly, $M$ in Eq. (94) is given, for large $\mu$, by

$$ M \approx \frac{1}{2}(\mu + l + \frac{d}{2}) \approx \frac{1}{2}(\mu + l). \qquad (105) $$

The subsequent analysis is different for $d = 1$ and $d \ge 2$ and we will analyse the two cases separately.

**Case** $d = 1$. In this case, there are only two allowed values of $l, l = 0$ and $l = 1$. Hence from Eqs. (104) and (105), it follows that, for large $\mu$, $m \approx M \approx \mu/2$ for both $l = 0$ and $l = 1$. Consequently from Eq. (101), one gets

$$ A_{m,M} \approx 2\mu \quad \text{and} \quad \sigma_{m,M} \approx 2\mu^{1/3}. \qquad (106) $$

Using $r_{max,l} = \sqrt{\lambda_{max}}$ from Eq. (96), it follows from Eq. (100) that, both for $l = 0$ and $l = 1$

$$ r_{max,l} = (A_{m,M} + \sigma_{m,M} \chi_2)^{1/2} \approx \sqrt{2\mu} + \frac{1}{\sqrt{2}} \mu^{-1/6} \chi_2, \qquad (107) $$

where we have expanded the square root for large $\mu$ and kept the two leading order terms. For $l = 0$ and $l = 1$, the CDF of $r_{max,l}$ thus has the scaling form

$$ P_l(w, m) = \text{Prob.}(r_{\text{max},l} \le w) = \mathcal{F}_2 \left( \sqrt{2} \mu^{1/6} (w - \sqrt{2}\mu) \right). \qquad (108) $$

Using this result in Eq. (99), it follows that in $d = 1$, the cumulative distribution of $r_{max}$ is the square of the Tracy-Widom distribution

$$ P(w, N) \approx \left[ \mathcal{F}_2 \left( \sqrt{2} \mu^{1/6} (w - \sqrt{2}\mu) \right) \right]^2, \qquad (109) $$
in the large N limit. This result has a natural interpretation which follows from Eq. (98). Indeed, Eq. (98) can be re-expressed in terms of the joint distribution of $x_{\text{min}} = \min\{x_1, \dots, x_N\}$ and $x_{\text{max}} = \max\{x_1, \dots, x_N\}$ as

$$P(w, N) = \text{Prob.}(x_{\text{max}} \le w, x_{\text{min}} \ge -w). \qquad (110)$$

While, for finite N, $x_{\text{max}}$ and $x_{\text{min}}$ are strongly correlated, they become independent for large N, leading to a decoupling of Eq. (110) into a product

$$P(w, N) \approx \text{Prob.}(x_{\text{max}} \le w) \text{Prob.}(x_{\text{min}} \ge -w). \qquad (111)$$

Since each of these probabilities is given, for large N, by the Tracy-Widom limiting distribution, one gets the square of the Tracy-Widom CDF in Eq. (109). Note that, even in the large N limit, this factorization does not hold for all w. It holds only when both the maximum and the minimum are close to their respective typical values, i.e. when $|w - \sqrt{2\mu}| \sim \mathcal{O}(\mu^{-1/6})$. In contrast when $|w - \sqrt{2\mu}| \sim \mathcal{O}(\sqrt{\mu})$, i.e., when both the maximum and the minimum are close to the center of the trap, this factorization no longer holds anymore. Instead $P(w, N)$ is described by a large deviation form that was computed exactly in [35].

**Case d ≥ 2.** In this case, the possible values of $l$ run from 0 to $\mu$, for large $\mu$. Substituting the large $\mu$ expansions given in Eqs. (104) and (105) in Eq. (101), expanding for large $\mu$ and for the moment keeping $l/\mu = \ell$ fixed, we obtain

$$A_{m,M} = \mu(1 + \sqrt{1 - \ell^2}) \quad \text{and} \quad \sigma_{m,M} = (2\mu)^{1/3} \frac{(1 + \sqrt{1 - \ell^2})^{2/3}}{(1 - \ell^2)^{1/6}} \quad \text{where} \quad \ell = \frac{l}{\mu}. \qquad (112)$$

Using $r_{\text{max},l} = \sqrt{\lambda_{\text{max}}}$ from Eq. (96), it follows from Eq. (100) that

$$r_{\text{max},l} = (A_{m,M} + \sigma_{m,M} \chi_2)^{1/2} \approx \sqrt{A_{m,M}} + \frac{1}{2} \frac{\sigma_{m,M}}{\sqrt{A_{m,M}}} \chi_2, \qquad (113)$$

where we have expanded the square root for large $\mu$ and kept the two leading order terms. This then gives

$$r_{\text{max},l} \approx C(\ell) \sqrt{\mu} + D(\ell) \mu^{-1/6} \chi_2, \qquad (114)$$

where

$$C(\ell) = (1 + \sqrt{1 - \ell^2})^{1/2} \quad \text{and} \quad D(\ell) = 2^{-2/3} \left[ \frac{1 + \sqrt{1 - \ell^2}}{1 - \ell^2} \right]^{1/6} \qquad (115)$$

where we recall that $\ell = l/\mu$. The CDF of $r_{\text{max},l}$ has thus the scaling form

$$P_l(w, m) = \text{Prob.}(r_{\text{max},l} \le w) = \mathcal{F}_2(s_\ell) \quad , \quad s_\ell = \frac{w - C(\ell) \sqrt{\mu}}{D(\ell) \mu^{-1/6}}, \qquad (116)$$

in the region of the typical fluctuations, where $s_\ell = \mathcal{O}(1)$ and we recall that $\mathcal{F}_2(s)$ is the CDF of the GUE TW distribution. We substitute this result (116) into the product in Eq. (82) and obtain

$$P(w, N) = \exp[S(w)] \quad \text{where} \quad S(w) = \sum_{l \ge 0} g_d(l) \ln \mathcal{F}_2 \left( \frac{w - C(l) \sqrt{\mu}}{D(l) \mu^{-1/6}} \right). \qquad (117)$$

For this formula to hold we need that the minimum value $s_{\text{min}} = \min_{\ell \in [0,1]} s_\ell$, of the argument $s_\ell$ of the TW distribution is in the range $s_{\text{min}} = \mathcal{O}(1)$ (see discussion below).

For $d \ge 2$, the sum over $l$ in Eq. (117) runs up to $\mu$ where $\mu$ is large. Hence the asymptotic large N analysis of $P(w, N)$ in this case is rather different from the $d=1$ case discussed above. It turns out that the leading contribution to the sum comes from the large $l$ regime (where $l = \mathcal{O}(\mu^{2/3})$ for large $\mu$, a fact to be verified a posteriori), where we can use the asymptotic form of $g_d(l)$ for large $l$ in Eq. (71). Furthermore, recalling that $l = \ell \mu$, it is clear that in the large $\mu$ limit, one can replace the sum over $l$ by an integral over $\ell$. This gives

$$S(w) \approx \frac{2 \mu^{d-1}}{\Gamma(d-1)} \int_0^1 d\ell \ell^{d-2} \ln \mathcal{F}_2 \left( \frac{w - C(\ell) \sqrt{\mu}}{D(\ell) \mu^{-1/6}} \right) \qquad (118)$$

where $C(\ell)$ and $D(\ell)$ are given in Eq. (115).

We must now analyze the argument $s_\ell = \mu^{2/3}(\hat{w} - C(\ell))/D(\ell)$ of the Tracy-Widom CDF, where we have defined $\hat{w} = w/\sqrt{\mu}$. Since $-\ln \mathcal{F}_2(s_\ell)$ is a positive decreasing function, the integral will be a priori dominated by the
minimum value of $s_\ell$, i.e. $s_{\min}$, as $\ell$ varies. The analysis depends on the value of $\hat{w}$, which we suppose fixed as $\mu$ becomes large. There are only two possible local minima for $s_\ell$ on the interval $\ell \in [0, 1]$, one at $\ell = 0$, where locally $s \simeq (2\mu)^{2/3}(2^{-1/6}(\hat{w} - \sqrt{2}) - \frac{1}{2^{1/68}}(\hat{w} - 2\sqrt{2})\ell^2)$ and one at $\ell = 1$ where $s \simeq (2\mu)^{2/3}(\hat{w} - 1)(1 - \ell^2)^{1/6}$. For $\hat{w} < 1$ there is only the minimum at $\ell = 0$. For $1 \le \hat{w} \le 2\sqrt{2}$ both minima exist, and the absolute minimum is at $\ell = 0$ for $\hat{w} < \sqrt{2}$, while it is at $\ell = 1$ for $\sqrt{2} < \hat{w} < 2\sqrt{2}$.

We now focus on the regime near the edge of the Fermi gas, $\hat{w} = w/\sqrt{\mu} = \sqrt{2}$, which turns out to contain the typical fluctuation regime (this can be anticipated from the results of section II B). We will see that the leading contribution comes from the small $\ell$ regime (of order $O(\mu^{-1/3})$). Let us examine the contribution of the minimum of $s_\ell$ at $\ell = 0$ (from the discussion of the previous paragraph). One can write the argument of $\mathcal{F}_2$ in Eq. (118) as

$$s_{\ell} = \frac{w - C(\ell) \sqrt{\mu}}{D(\ell) \mu^{-1/6}} \approx \tilde{w} + \frac{1}{4} \ell^2 \mu^{2/3} + \mathcal{O}(\ell^4) \quad \text{where} \quad \tilde{w} = (w - \sqrt{2\mu}) \sqrt{2} \mu^{1/6}. \qquad (119)$$

Clearly one must rescale $\ell = 2v\mu^{-1/3}$ in the integral over $\ell$ and one obtains in the large $\mu$ limit

$$S(w) \approx \frac{2^d \mu^{\frac{2(d-1)}{3}}}{\Gamma(d-1)} \int_0^{+\infty} dv \, v^{d-2} \ln \mathcal{F}_2(\tilde{w} + v^2) \qquad (120)$$

This result is correct as long as $\hat{w} - \sqrt{2}$ vanishes at large $\mu$, i.e. $|\tilde{w}| \ll \mu^{2/3}$. For other values of $\hat{w}$ one needs to use, instead of the typical fluctuation formula given in Eq. (116), a large deviation formula, which goes beyond the scope of this paper.

Let us now study the regime of typical fluctuations. Denoting, as in Eqs. (43)-(45), $\tilde{w} = a_\mu + z/(2\sqrt{a_\mu})$ and anticipating that $a_\mu$ is large and $z = O(1)$, we can now use the asymptotic tail of $\mathcal{F}_2(s)$ for large $s$ in Eq. (103), inserting $s = a_\mu + z/(2\sqrt{a_\mu}) + v^2$ in (120), one finds for large $\mu$

$$S(w) \approx - \frac{2^d \mu^{\frac{2(d-1)}{3}}}{16\pi\Gamma(d-1) a_\mu^{3/2}} \int_0^{+\infty} dv \, v^{d-2} e^{-\frac{4}{3}a_\mu^{3/2} - 2\sqrt{a_\mu}v^2 - z}. \qquad (121)$$

We can now substitute, as in (43)-(45), $\tilde{w} = a_\mu + z/(2\sqrt{a_\mu})$ where $a_\mu$ is large and $z = O(1)$. Expanding, the leading behavior at large $\mu$ of the expression (121) leads to evaluate the integral

$$\int_0^{+\infty} dv \, v^{d-2} e^{-2\sqrt{a_\mu}v^2} = 2^{-\frac{d+1}{2}} a_\mu^{\frac{1-d}{4}} \Gamma\left(\frac{d-1}{2}\right), \qquad (122)$$

leading exactly to

$$S(w) \approx -B_d e^{-\frac{4}{3} a_\mu^{3/2} + \frac{2(d-1)}{3} \ln \mu - \frac{(d+5)}{4} \ln(a_\mu) - 2 \sqrt{a_\mu} b_\mu z}, \qquad (123)$$

where $B_d$ was given in (42). Hence, we arrive exactly at the result in Eq. (44), and (49) in section II B, i.e. a Gumbel distribution for $r_{\max}$, although by a quite different method. The present calculation further shows that the angular momentum which contributes to the value of $r_{\max}$ are of order $l = \ell\mu \sim \mu^{2/3}$ (see the interpretation below).

Let us comment on the contribution of the second minimum of $s_\ell$ at $\ell = 1$, which corresponds to high angular momentum states, and turns out to be negligible for $\hat{w} = \sqrt{2}$. A naive analysis of the continuum integral (118) would suggest an additional contribution to $S(w)$ corresponding to values of $\ell = 1 - O(\mu^{-4})$. However in the original sum (117), $l$ is an integer which runs up to $l_{\max} = \text{int}(\mu - d/2)$. Hence the continuum integral must be cutoff at values of $\ell = 1 - O(\mu^{-1})$, which corresponds to $s_{\min} \sim \sqrt{\mu}$, hence leading to an exponentially small contribution to $S(w)$, i.e. $\sim O(\mu^{d-2}e^{-\mu^{3/4}}) \ll \mu^{\frac{2(d-1)}{3}}$. In fact one can directly analyze the contribution coming from $l = l_{\max}$, corresponding to states with $n = 0$ and having a degeneracy $\sim \mu^{d-2}$ (see Fig. 2). The associated radial wavefunction $\chi_{0,l_{\max}}(r)$ given by (61) can be shown to be peaked around $r^* \simeq \sqrt{\mu} \ll$ which is far from the trap center but is also far from the edge of the Fermi gas (recall $r_{\text{edge}} = \sqrt{2\mu}$). One can calculate explicitly for this state, using formula (84) (which is a simple one-dimensional determinant as $m_l = 1$) the probability

$$P_{l_{\max}}(w, m_{l_{\max}} = 1) = 1 - \frac{\int_{w}^{+\infty} dr \, r^{2l_{\max}+d-1} e^{-r^2}}{\int_{0}^{+\infty} dr \, r^{2l_{\max}+d-1} e^{-r^2}} \approx 1 - O(e^{-2\mu}), \qquad (124)$$

leading to an exponentially small contribution to $S(w)$, i.e. $\sim O(\mu^{d-2}e^{-2\mu}) \ll \mu^{\frac{2(d-1)}{3}}$. Hence there is no additional contribution due to these maximal angular momentum states, in agreement with the analysis above.
# IV. EXTENSIONS TO OTHER OBSERVABLES AND FINITE TEMPERATURE

In this section we give some further probabilistic interpretations as well as new results. We study counting statistics near the edge, extreme value statistics in non-spherical domains, the distribution of $x_{max}$, and present some extensions to finite temperature.

## A. Poisson statistics near the expected maximal radial distance

Let us define the number of fermions $N_w$ outside the sphere centered at the origin and of radius $w$ in the ground state. From the determinantal structure of the quantum probability of $N$ fermions at zero temperature (19), we know that the Laplace transform of the distribution of $N_w$ can be written as a Fredholm determinant

$$ \langle e^{-pN_w} \rangle_0 = \mathrm{Det}(I - (1 - e^{-p})I_w K_\mu). \qquad (125) $$

We will restrict our analysis to values of $w$ near the edge $w \approx r_{edge}$ where, if $p = O(1)$ and in $d > 1$, we can use the single trace approximation as in Section II B. In addition we consider the region around the average value of $r_{max}$, hence we set

$$ w = \sqrt{2\mu} + \frac{1}{2\mu^{1/6}} \left( a_{\mu} + \frac{z}{\sqrt{2a_{\mu}}} \right), \qquad (126) $$

where $z$ is a parameter of order unity. Following the same steps as in that Section we easily arrive at

$$ \langle e^{-pN_w} \rangle_0 = \exp((e^{-p} - 1)e^{-z}). \qquad (127) $$

This means that $N_w$ has a Poissonian statistics with parameter $e^{-z}$, i.e.

$$ \mathrm{Prob}(N_w = k) = \frac{e^{-kz}}{k!} \exp(-e^{-z}), \qquad (128) $$

where the random variable is $N_w$ and $z$ a parameter. Since one has $\langle N_w \rangle_0 = e^{-z}$ the Poisson distribution crosses over to a Gaussian distribution for $z < 0$ and $|z| \gg 1$. This is the case in the region close to $r_{edge}$, where $\tilde{w} = O(1) \ll a_\mu$.

## B. Other domains

Studying systems in higher dimensions opens up the possibility to study a wider variety of extreme value problems than in one dimension. For instance one can ask what is the distribution of the farthest fermion from the origin when one looks within a certain solid angle, or what is the distribution of the maximal coordinate in a given direction. Here we address these issues.

Let us come back to the single trace approximation in Eq. (36) for the probability that $r_{max} < w$, namely

$$ \mathrm{Prob}(r_{\mathrm{max}} \le w) \approx \exp \left[ - \int_{\mathbf{x} \in \mathcal{D}} K_{\mu}(\mathbf{x}, \mathbf{x}) d\mathbf{x} \right] \approx \exp \left[ -S_d \left( \frac{r_{\mathrm{edge}}}{w_N} \right)^{d-1} \int_{\tilde{w}}^{\infty} F_d(s) ds \right], \qquad (129) $$

where here the domain $\mathcal{D}$ is the region $|\mathbf{x}| > w$, and we have chosen a $w$ such that $w - r_{edge} = O(w_N)$. One can interpret the second equality, roughly, as saying that there are a number of order $N_{eff} \sim S_d(r_{edge}/w_N)^{d-1}$ independent regions contributing to the maximal radius (naturally leading then to a Gumbel distribution). Said differently, it means that along the surface of the Fermi droplet, i.e. in the angular direction, the positions of the fermions within the domain $\mathcal{D}$, are correlated over distances of order $w_N$. This is only a rough interpretation, a more detailed and precise one was presented in the previous section, but it is sufficient for our purpose here.

It is then reasonable to assume that the same approximation should hold for more general domains, as long as the number of independent regions $N_{eff}$ is large as $N \to +\infty$, leading again to a Gumbel distribution (this condition also ensures that the higher power traces in the expansion of the Fredholm determinant (32) remain subdominant). Let us for instance consider a domain $\mathcal{D}$ defined as the region $|\mathbf{x}| > w(\boldsymbol{\theta})$ where $\boldsymbol{\theta}$ is an angular direction (i.e. a $(d-1)$-dimensional angular vector, e.g. a simple angle $\theta \in [0, 2\pi[$ for $d=2$). Here $w(\boldsymbol{\theta})$ is an arbitrary curve, but to be able to use the edge scaling form of the kernel we will choose it such as $w(\boldsymbol{\theta}) \ge r_{edge} + O(w_N)$ and we define

$$ \tilde{w}(\boldsymbol{\theta}) = \frac{w(\boldsymbol{\theta}) - r_{\mathrm{edge}}}{w_N} = (w(\boldsymbol{\theta}) - \sqrt{2\mu}) \sqrt{2\mu}^{1/6}. \qquad (130) $$
<figure><img src="image_4.png" /><figcaption>FIG. 3. Illustration of the domain $\mathcal{D}$, in $d = 2$, corresponding to $w(\theta) = w$ in the angular sector $\theta \in [0, \theta_m]$, and elsewhere $\tilde{w}(\theta) = +\infty$. The hole probability associated with $\mathcal{D}$ is the probability that $r_{\text{max}} \le w$ restricted to the fermions in the angular sector $[0, \theta_m]$.</figcaption></figure>

We then ask what is the probability that there are no fermions in the domain $\mathcal{D}$, i.e. the “hole probability” associated to the domain $\mathcal{D}$. One example in $d = 2$ is $\tilde{w}(\theta) = \tilde{w}$ in the angular sector $\theta \in [0, \theta_m]$, and else $\tilde{w}(\theta) = +\infty$, and this hole probability then becomes the probability that $r_{\text{max}} \le w$ restricted to the fermions in the angular sector $[0, \theta_m]$ (see Fig. 3). In general it can be expressed as a Fredholm determinant as in (30), (31), with $I_w(\mathbf{x})$ the indicator function of the domain $\mathcal{D}$. Keeping again the leading single trace term we find (defining $N_{\mathcal{D}}$ the number of fermions in the domain $\mathcal{D}$)

$$ \text{Prob}(N_{\mathcal{D}} = 0) \approx \exp \left[ - \left( \frac{r_{\text{edge}}}{w_N} \right)^{d-1} \int d\theta \int_{\tilde{w}(\theta)}^{\infty} F_d(s) ds \right], \quad (131) $$

and similar manipulations as in Section IIB lead to a Gumbel distribution for a suitably scaled $\tilde{w}$ (with different constants $a_\mu$ and $b_\mu$).

Let us give an example of application in $d = 2$ (another example being given in the next subsection). First consider the probability, denoted $P_{\theta_m}(w, N)$, that $r_{\text{max}} \le w$, restricted to the fermions in the angular sector $[0, \theta_m]$. The effective number of independent regions is $N_{\text{eff}} \sim \theta_m \mu^{2/3}$. Hence if the angular sector is sufficiently large, i.e. if $\theta_m \gg \mu^{-2/3} \sim N^{-1/3}$ one obtains a Gumbel distribution, $P_{\theta_m}(w, N) \approx \exp(-e^{-z})$ where $\tilde{w} = a_\mu + \frac{1}{2\sqrt{a_\mu u}}z$, with a parameter $a_\mu$ solution of

$$ \frac{2^{-7/2}}{\sqrt{\pi}} \theta_m e^{-\frac{4}{3} a_{\mu}^{3/2} + \frac{2}{3} \ln \mu - \frac{7}{4} \ln(a_{\mu})} = 1, \quad (132) $$

leading to $a_\mu \approx [\frac{1}{2} \ln(\mu \theta_m^{3/2})]^{\frac{2}{3}}$. If the angular sector becomes too small, i.e., $\theta_m$ of order, or smaller than $\mu^{-2/3} \sim N^{-1/3}$ we expect a crossover to a Tracy-Widom regime, which would be quite interesting to study, but is beyond the scope of the present paper. Note that the scale $\theta_m \sim \mu^{-2/3}$ is consistent with the angular momentum scale $l \sim \theta_m^{-1} \sim \mu^{2/3}$ which was shown in the previous Section to control the value of $r_{\text{max}}$.

Note finally that the arguments in this section have been mostly qualitative, and to make them rigorous (and study the crossover to the $d = 1$ limit) would require a general analysis of the higher traces for a larger class of domains $\mathcal{D}$ which goes beyond the scope of the present paper.

### C. Probability distribution of $x_{\text{max}}$ in $d = 2$

Another interesting question is to calculate the PDF of the first coordinate, i.e. its abscissa denoted by $x_{\text{max}}$, of the rightmost fermion in $d = 2$ dimensions. We start with a first method, using results of the previous subsection. This corresponds to the choice $w(\theta) = w/\cos(\theta)$. One can check a-posteriori that one can restrict the analysis to small
angles $\theta$, hence we write $w(\theta) = w(1 + \frac{\theta^2}{2})$. Thus we obtain

$$ \text{Prob.}(x_{\text{max}} \le w) \approx \exp \left[ -2\mu^{2/3} \int d\theta \int_{\tilde{w}+\mu^{2/3}\theta^2}^{\infty} F_d(s) ds \right], \quad (133) $$

where we have used the definition (130). Expanding as $\tilde{w} = a_\mu + \frac{1}{2\sqrt{a_\mu}}z$ we obtain again a Gumbel distribution for the variable $z$, i.e. $\text{Prob.}(x_{\text{max}} \le w) = e^{-e^{-z}}$, with now the following condition to determine $a_\mu$

$$ \tilde{B}_2(\mu) e^{-\frac{4}{3} a_{\mu}^{3/2} + \frac{2}{3} \ln \mu - \frac{7}{4} \ln(a_{\mu})} = 1, \quad \tilde{B}_2(\mu) = \frac{1}{(8\pi)^{3/2}} \int_{-\infty}^{\infty} d\theta e^{-2a_{\mu}^{1/2} \mu^{2/3} \theta^2} = \frac{1}{32\pi} a_{\mu}^{-1/4} \mu^{-1/3}, \quad (134) $$

which leads to $-\frac{4}{3} a_\mu^{3/2} + \frac{1}{3} \ln \mu - 2 \ln a_\mu - \ln(32\pi) = 0$. Hence to leading order one has

$$ a_{\mu} \approx \left( \frac{1}{4} \ln \mu \right)^{2/3}. \quad (135) $$

Note that the angles involved in this maximization are $\theta \sim \theta_\mu = a_\mu^{-1/4} \mu^{-1/3}$ which justifies the expansion of the $\cos(\theta)$ for small $\theta$ (i.e., its approximation by a parabola). The maximum $x_{\text{max}}$ is thus attained for fermions which deviate at most of order $\theta_\mu$ from the axis (up to powers of $\ln \mu$). Also the number of independent regions is $N_{\text{eff}} = \mu^{1/3} a_\mu^{-1/4} \gg 1$ as $N$ becomes large, which a posteriori justifies the calculation. To summarize, we find

$$ x_{\text{max}} = \sqrt{2\mu} + \frac{1}{\mu^{1/6}} \left[ \frac{1}{\sqrt{2}} a_{\mu} + \frac{1}{2\sqrt{2} a_{\mu}} Z_G \right], \quad (136) $$

where $\text{Prob.}[Z_G \le z] = \exp[-e^{-z}]$ and $a_\mu$ is given by (135). Comparing with the similar formula for $r_{\text{max}}$, in Eq. (49), (50) we see that $x_{\text{max}}$ is smaller than $r_{\text{max}}$ but its fluctuations are larger.

We will now show that this result can be obtained from an exact calculation, similar to the one of Section III, which we present here only in $d=2$ for simplicity – the generalization to higher $d$ being straightforward. To study $x_{\text{max}}$ it is more convenient to use the eigenbasis of the single particle problem

$$ \phi_{n_x, n_y}(x, y) = \phi_{n_x}(x)\phi_{n_y}(y) \quad , \quad \phi_k(x) = \left( \frac{1}{\sqrt{\pi 2^k k!}} \right)^{1/2} H_k(x)e^{-x^2/2}, \quad (137) $$

where $\phi_k(x)$ are the eigenstates of the $1d$ harmonic oscillator, and $H_k$ the Hermite polynomial of degree $k$. As in formula (80) we have

$$ \begin{aligned} \text{Prob.}(x_{\text{max}} \le w) &= \det_{1 \le i,j \le N} \left[ \delta_{ij} - \int_{x \ge w} dx dy \phi_{n_{x,i},n_{y,i}}(x,y) \phi_{n_{x,j},n_{y,j}}(x,y) \right] \\ &= \det_{1 \le i,j \le N} \left[ \delta_{ij} - \delta_{n_{y,i},n_{y,j}} \int_{x \ge w} dx \phi_{n_{x,i}}(x) \phi_{n_{x,j}}(x) \right] \end{aligned} \quad (138) $$

where the $(n_{x,i}, n_{y,i})$, $i = 1, \dots, N$ denote the occupied states (note that in Eq. (138) we have used that the eigenfunctions $\phi_k(x)$ are real). In the $N$-body ground state all single particle states with energy $\epsilon(n_x, n_y) = n_x + n_y + 1 \le \mu$ are occupied, leading to the relation between $N$ and the Fermi energy $\mu$: $N = \frac{1}{2}\mu(\mu + 1)$ (we consider $\mu$ to be equal to the energy of the last occupied state, hence it is an integer).

As in the previous calculation (see section III) this determinant (138) has a block structure and hence one can write

$$ \text{Prob.}(x_{\text{max}} \le w) = \prod_{n_y \ge 0}^{\mu-1} [P_{n_y}(w, \mu - 1 - n_y)] \quad (139) $$

where here

$$ P_{n_y}(w, N_{n_y}) = \det_{1 \le i,j \le N_{n_y}} \left[ \delta_{ij} - \int_{x \ge w} dx \phi_{n_{x,i}}(x) \phi_{n_{x,j}}(x) \right] \quad (140) $$

is the cumulative distribution function (CDF) of the position of the right-most fermion among $N_{n_y}$ fermions in a one-dimensional harmonic potential. This CDF has been studied in Ref. [8] – as well as in the RMT literature [13] – and is known to take the form of a Tracy-Widom distribution at large $N_y$, i.e. one has

$$ x_{\text{max},n_y} = \sqrt{2(\mu - 1 - n_y)} + \frac{1}{\sqrt{2}}(\mu - 1 - n_y)^{-1/6} \chi_2, \quad (141) $$
<figure><img src="image_5.png" /><figcaption>FIG. 4. Illustration of the position $x_{max}$ of the rightmost fermion for a two-dimensional system confined in a harmonic trap, as studied in [see Eq. (139)]. In the large $N$ limit, the PDF of $x_{max}$ properly centered and scaled converges also to a Gumbel distribution.</figcaption></figure>

where the CDF of $\chi_2$ is $\mathcal{F}_2(s)$ (to be compared with formula (113) for the study of $r_{max}$). Setting $n_y = \ell\mu$ in the limit of large $\mu$ and following similar steps as in Section III, one finds $\text{Prob}(x_{max} \le w) \simeq \exp(S(w))$ where $S(w)$ can be written as a sum over the states

$$S(w) = \mu \int_{0}^{1} d\ell \ln \mathcal{F}_{2} \left( (w - \sqrt{2\mu(1-\ell)}) \sqrt{2\mu}^{1/6} (1-\ell)^{1/6} \right), \quad (142)$$

where we have made the approximation $\sum_{n_y=0}^{\mu-1} \simeq \mu \int_0^1 d\ell$. Inserting the right tail of the Tracy-Widom distribution, we find again that the region which dominates the integral is for small $\ell$. Hence, expanding in $\ell$ close to $\ell = 0$, the argument of $\mathcal{F}_2$ can be replaced by $\tilde{w} + \ell\mu^{2/3}$. A similar calculation to that in Section III then allows us to recover exactly the result of Eq. (136), with exactly the same condition determining $a_\mu$, leading to the estimate in (135).

### D. Higher traces and the regime $\tilde{w} = \mathcal{O}(1)$

It is also interesting to ask what is the effect of the higher order traces in (34). As we have hinted at above they only contribute to subdominant corrections in the determination of the distribution of the typical fluctuations of $r_{max}$ [see Eq. (49)], when $\tilde{w} \sim a_\mu \gg 1$. However these terms become important if one studies the regime $\tilde{w} = \mathcal{O}(1)$, as we now discuss.

In Ref. [8, 11] we have shown that for $\mathbf{x}, \mathbf{y}$ near a point $r_{edge}\mathbf{n}$, where $\mathbf{n}$ is a unit vector (on the unit sphere), the kernel takes the scaling form

$$K_{\mu}(\mathbf{x}, \mathbf{y}) \approx \frac{1}{w_N^d} K_{d, \mathbf{n}}^{\text{edge}} \left( \frac{\mathbf{x} - r_{\text{edge}}\mathbf{n}}{w_N}, \frac{\mathbf{y} - r_{\text{edge}}\mathbf{n}}{w_N} \right) \quad (143)$$

where the $d$-dimensional edge kernel can be written as

$$K_{d, \mathbf{n}}^{\text{edge}}(\mathbf{a}, \mathbf{b}) = \int \frac{d\mathbf{q}_t}{(2\pi)^{d-1}} e^{-i\mathbf{q}_t \cdot (\mathbf{a}_t - \mathbf{b}_t)} K_{\text{Ai}}(a_n + \mathbf{q}_t^2, b_n + \mathbf{q}_t^2), \quad (144)$$

in terms of the standard Airy kernel, $K_{\text{Ai}}(a, b) = \int_0^{+\infty} du \text{Ai}(a+u)\text{Ai}(b+u)$. In Eq. (144) and below we decompose any vector $\mathbf{a} = \mathbf{a}_t + a_n\mathbf{n}$ into a normal (scalar) component $a_n = \mathbf{a} \cdot \mathbf{n}$ and a tangential component $\mathbf{a}_t$ (a ($d-1$)-dimensional vector). We have indicated by an additional subscript that the edge kernel explicitly depends on $\mathbf{n}$. Let us note the following interesting property

$$\int d\mathbf{b}_t K_{d, \mathbf{n}}^{\text{edge}}(\mathbf{a}, \mathbf{b}) K_{d, \mathbf{n}}^{\text{edge}}(\mathbf{b}, \mathbf{c}) = \int \frac{d\mathbf{q}_t}{(2\pi)^{d-1}} e^{-i\mathbf{q}_t \cdot (\mathbf{a}_t - \mathbf{c}_t)} K_{\text{Ai}}(a_n + \mathbf{q}_t^2, b_n + \mathbf{q}_t^2) K_{\text{Ai}}(b_n + \mathbf{q}_t^2, c_n + \mathbf{q}_t^2) \quad (145)$$
which allows (by iteration, see below and Appendix A) to calculate any power of the kernel in terms of Airy functions.

Using this scaling form (143) one can reduce the problem of calculating the trace in Eq. (35) to free integrations over $\mathbf{a}_2, \dots, \mathbf{a}_n$ together with an integral over $\mathbf{n}$ (over the unit sphere) and over $a_{1,n}$ the radial part of $\mathbf{a}_1$. Hence it is easy to see that (see Appendix A for details)

$$ P(w, N) \approx \exp \left( - \left( \frac{r_{\text{edge}}}{w_N} \right)^{d-1} \sum_{p=1}^{\infty} \frac{1}{p} \int d\mathbf{n} \operatorname{Tr} \left[ \delta_{\mathbf{a}_{1,t}, \mathbf{0}} (\tilde{I}_{\tilde{w}, \mathbf{n}} K_{d, \mathbf{n}}^{\text{edge}})^p \right] \right) \quad (146) $$

$$ = \exp \left( \left( \frac{r_{\text{edge}}}{w_N} \right)^{d-1} \int d\mathbf{n} \ln \operatorname{Det}(I - \delta_{\mathbf{a}_{1,t}, \mathbf{0}} \tilde{I}_{\tilde{w}, \mathbf{n}} K_{d, \mathbf{n}}^{\text{edge}}) \right), \quad (147) $$

where $\int d\mathbf{n}$ denotes the integral over the unit sphere, $\delta_{\mathbf{a}_{1,t}, \mathbf{0}}$ denotes the projector on the subspace $\mathbf{a}_{1,t} = \mathbf{0}$, and the trace is over all vectors $\mathbf{a}_1, \dots, \mathbf{a}_n$. In the last line the traces have been summed into a Fredholm determinant. If one studies $r_{\text{max}}$, then the operator $\tilde{I}_{\tilde{w}, \mathbf{n}}(\mathbf{a})$ is the projector on $a_n \ge \tilde{w}$, but the above expression is valid for more general domains (such as studied in Section IV B). In the case of the calculation of $r_{\text{max}}$, using rotational invariance it is easy to see (see Appendix A) that all the traces in (146) are independent of $\mathbf{n}$ hence one can replace $\int d\mathbf{n}$ by $S_d$ and pick one particular fixed direction $\mathbf{n}$. The calculation of the traces then uses an iteration of (145)

$$ \operatorname{Tr} \left[ \delta_{\mathbf{a}_{1,t}, \mathbf{0}} (\tilde{I}_{\tilde{w}, \mathbf{n}} K_{d, \mathbf{n}}^{\text{edge}})^p \right] = \int \frac{d\mathbf{q}_t}{(2\pi)^{d-1}} \operatorname{Tr} \left[ (P_{\tilde{w}+\mathbf{q}_t^2} K_{Ai})^p \right], \quad (148) $$

where $P_s$ is the standard notation for the projector on $[s, +\infty[$. Summing up the traces we obtain the formula for the probability that $r_{\text{max}} \le w$:

$$ P(w, N) \approx \exp \left( S_d \left( \frac{r_{\text{edge}}}{w_N} \right)^{d-1} \int \frac{d\mathbf{q}_t}{(2\pi)^{d-1}} \ln \operatorname{Det}(I - P_{\tilde{w}+\mathbf{q}_t^2} K_{Ai}) \right) \quad (149) $$

$$ = \exp \left( S_d \left( \frac{r_{\text{edge}}}{w_N} \right)^{d-1} \int \frac{d\mathbf{q}_t}{(2\pi)^{d-1}} \ln \mathcal{F}_2(\tilde{w} + \mathbf{q}_t^2) \right), \quad (150) $$

where $\mathcal{F}_2(s) = \operatorname{Det}(I - P_s K_{Ai})$ is the GUE Tracy-Widom CDF. This expression is valid for $N$ large, and for $\tilde{w} = \mathcal{O}(1)$, as well as in the regime $\tilde{w} \sim a_\mu \gg 1$ studied previously. Hence this expression shows the large deviation form for the PDF of $r_{\text{max}}$, when $\tilde{w} = \mathcal{O}(1)$, which is an intermediate left deviation regime for the present problem (note that there exists another left large deviation regime, deeper inside the Fermi cloud, whose study is left for the future). The exponent in (150) is the associated rate function. Here we have obtained it exactly in terms of the GUE-TW function $\mathcal{F}_2(s)$, in the case of a spherically symmetric domain.

Performing the angular integral in (150), and denoting $|\mathbf{q}_t| = v$ it is immediate to see, using the identities $S_d S_{d-1}/(\pi)^{d-1} = 2^d/\Gamma(d-1)$ and $r_{\text{edge}}/w_N = 2\mu^{2/3}$, that the formula (150) is identical to Eq. (120), obtained in Section III D by quite different methods. There it was studied only in the typical fluctuation regime $\tilde{w} \simeq a_\mu \gg 1$, but the present analysis shows that these formula holds down to the regime $\tilde{w} = \mathcal{O}(1)$.

The explicit calculation of the rate function for more general domains in this regime $\tilde{w} = \mathcal{O}(1)$, i.e. of the reduced Fredholm determinant for non-spherically symmetric domains, remains challenging. The analysis of the previous section shows, that a simplification should also arise in the calculation of $x_{\text{max}}$. Finally, in the region $\tilde{w} \sim a_\mu$ this Fredholm determinant can be approximated by the leading single trace, and one recovers the typical fluctuations studied in this paper (see Appendix A for further details).

## E. Finite temperature

We now briefly sketch the effect of a finite temperature on the distribution of the maximal radius, denoted by $r_{\text{max}}(T)$. As above, we restrict ourselves to the case of the harmonic oscillator. As was shown in [11] the corresponding hole probability can be calculated at finite temperature. For this purpose, one uses the equivalence of the canonical (fixed $N$) and grand-canonical (fixed chemical potential $\tilde{\mu}$) ensembles in the large $N$ limit, and the fact that in the grand-canonical ensemble the JPDF of the positions of the fermions form a determinantal process [11, 36]. As a result the formula (34) is still valid but $K_\mu$ is replaced by the finite temperature kernel, described in [11]. Since here we study $r_{\text{max}}(T)$, i.e., the edge of the Fermi gas, we can again use the scaling form of the finite temperature kernel near the edge, which is valid in the regime where the temperature $T \sim N^{1/(3d)}$: this is the regime on which we focus from now on. More precisely one defines the dimensionless parameter [9, 11]

$$ b = (\Gamma(1+d)N)^{1/(3d)}/T. \quad (151) $$
We now use the same (single trace) approximation as we used at $T = 0$ (whose validity is discussed below) and we obtain, following similar steps to those in Section II B

$$P(w, N) \approx \exp \left[ -S_d \left( \frac{r_{\text{edge}}}{w_N} \right)^{d-1} \int_{\tilde{w}}^{\infty} F_{d,b}^{\text{edge}}(z) dz \right]. \quad (152)$$

As before, we have written again $\tilde{w} = \frac{w-r_{\text{edge}}}{w_N} = (w - \sqrt{2\mu}) \sqrt{2} \mu^{1/6}$ where here and below $\mu$ is the Fermi energy, i.e. the $T = 0$ chemical potential, as defined in Section II [see Eq. (14)]. This formula involves the finite temperature version of the scaling function of the density, obtained in [11] as

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

where

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

is the $d=1$ finite temperature density scaling function. Clearly we again need the large positive $\tilde{w}$ asymptotics of Eq. (152), which requires the large argument asymptotics of the function $F_{d,b}(z)$ (right tail), which was not given in [11]. Let us recall that the right tail of $F_{1,b}(z)$ was obtained in [11, 37] and shown to depend on the parameter $\tilde{s} = s/b^2$ as follows

$$F_{1,b}^{\text{edge}}(s) \approx \begin{cases} \frac{1}{4b^2 \sqrt{\tilde{s}} \sin(2\pi\sqrt{\tilde{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 (155)$$

Hence it exhibits a transition between a stretched exponential tail to a pure exponential decay for $s > s_c = b^2/4$. Thus, for a fixed value of the reduced temperature $b$ (not necessarily large), the decay is always exponential. The pre-exponential factor further exhibits a crossover from the two limiting cases indicated above, in the vicinity of $\tilde{s} = 1/4$ as given in [37]. Note that the first regime in Eq. (155) for $\tilde{s} \to 0$, when inserted in (153), gives back the $T=0$ behavior (40).

Let us start with the “thermal” regime where we can use the second form in (155), i.e. $s > \frac{b^2}{4}$. Inserting into (153) we find the right tail

$$F_{d,b}^{\text{edge}}(s) \approx \frac{1}{(4\pi b)^{d/2}} \exp\left(-bs + \frac{b^3}{12}\right), \quad (156)$$

which gives

$$P(w, N) \approx \exp \left[ -S_d \left( \frac{r_{\text{edge}}}{w_N} \right)^{d-1} \frac{1}{(4\pi)^{d/2} b^{1+\frac{d}{2}}} e^{-b\tilde{w}+\frac{b^3}{12}} \right], \quad (157)$$

where we recall that $\frac{r_{\text{edge}}}{w_N} = 2\mu^{2/3} = 2(\Gamma(d+1)N)^{2/(3d)}$. Hence we obtain again a Gumbel distribution in terms of a variable $z$ now defined as

$$P(w, N) \approx e^{-e^{-z}}, \quad \tilde{w} = a_{\mu} + \frac{1}{b}z \quad (158)$$

where $a_\mu$ is here determined by the condition

$$A_d e^{-ba_\mu + \frac{2}{3}(d-1)\ln\mu + \frac{b^3}{12} - (1+\frac{d}{2})\ln b} = 1, \quad (159)$$

with $A_d = \frac{S_d 2^{d-1}}{(4\pi)^{d/2}} = 1/\Gamma(d/2)$. This condition (159) yields the leading order behavior

$$a_{\mu} = a_{\mu,b} \approx \frac{2}{3b}(d-1)\ln\mu + \frac{b^2}{12} - \frac{d+2}{2b}\ln b + \frac{1}{b}\ln A_d, \quad (160)$$
which is our final result for the distribution of $r_{\max}(T)$ in the thermal regime, i.e., of fixed $b$ and large $N$ (i.e. large $\mu$)

$$r_{\max}(T) = \sqrt{2\mu} + \frac{1}{\sqrt{2\mu^{1/6}}} \left( a_{\mu,b} + \frac{1}{b} Z_G \right) \qquad (161)$$

where Prob.[$Z_G \le z$] = $\exp[-e^{-z}]$. We note that this result remains valid for $d=1$ in the limit $b \ll 1$, and exactly matches with the result obtained in [11] (discarding the term $b^2/12$ which is clearly subdominant in this limit). To check this we used the fact – see (111) – that in $d=1$, the PDF of $r_{\max}(T)$ is given by the maximum of two independent Gumbel random variables, each obtained in Eq. (156) of [11].

Clearly, comparing this result in Eqs. (160) and (161) with the result obtained at $T=0$ and given in Eqs. (49) and (50), we see that there is a crossover at very low temperature $b = b^*$ from the thermal regime to the $T=0$ regime, with

$$b^* \sim \frac{4}{3} \left( \frac{d-1}{2} \ln \mu \right)^{1/3} . \qquad (162)$$

This crossover occurs between two Gumbel distributions with different parameters. The study of the details of this crossover is left for the future.

Let us also indicate the result for the distribution of $x_{\max}(T)$ at finite temperature in $d=2$. By a similar calculation, using similar arguments as in Section IV C, we find

$$x_{\max}(T) = \sqrt{2\mu} + \frac{1}{\sqrt{2\mu^{1/6}}} \left( a_{\mu,b} + \frac{1}{b} Z_G \right) \qquad (163)$$

where Prob.[$Z_G \le z$] = $\exp[-e^{-z}]$ is a Gumbel variable, with

$$a_{\mu} = a_{\mu,b} \approx \frac{1}{3b} \ln \mu + \frac{b^2}{12} - \frac{5}{2b} \ln b - \frac{1}{b} \ln(2\sqrt{\pi}) . \qquad (164)$$

Note that we find that the amplitude of the fluctuations of $x_{\max}(T)$ and $r_{\max}(T)$ are identical but that on average $r_{\max}(T) > x_{\max}(T)$, as expected.

Finally, recalling from [11] (section VII-D) that the formula (144) extends to any value of $b$, one easily sees that the expression in (150) extends to finite temperature, at fixed $b$, as follows

$$P(w, N) \approx \exp \left( S_d \left( \frac{r_{\text{edge}}}{w_N} \right)^{d-1} \int \frac{d\mathbf{q}_t}{(2\pi)^{d-1}} \ln \mathcal{F}_{2,b}(\tilde{w} + \mathbf{q}_t^2) \right), \qquad (165)$$

where here we denote $\mathcal{F}_{2,b}(s) = \text{Det}(I - P_s K_b^{\text{edge}})$ the finite temperature generalization of the GUE Tracy-Widom CDF, where $K_b^{\text{edge}}$ is given in section V-C of [11].

## V. DISCUSSION

In this paper, we have investigated analytically the distribution of the maximal radial distance from the trap center of $N$ non-interacting fermions in a $d$-dimensional harmonic trap at $T=0$. We have shown that, in the large $N$ limit, the cumulative distribution, appropriately centered and scaled, converges to the Gumbel distribution for $d > 1$. In $d=1$, we have shown that the limiting cumulative distribution is given by the square of the Tracy-Widom distribution for GUE. If one thinks of the dimension $d$ as a continuous parameter, clearly the distribution crosses over from the squared Tracy-Widom to the Gumbel form. It would be interesting to investigate the precise crossover function that interpolates between these two distinct forms, as $d$ crosses the critical value $d=1$.

How universal is the large $N$ scaling form of the cumulative distribution with respect to the shape of the trapping potential? We have shown that the limiting form of the cumulative distribution can be derived solely from the asymptotic properties of the edge kernel. However it has been recently shown that the edge kernel, properly centered and scaled, is universal, i.e., independent of the details of the trapping potential, for a broad class of confining potentials. Consequently, the limiting distribution of the maximal radial distance is also universal, with respect to this class of trapping potentials. Note that a different extreme value observable is the smallest distance to the center of the trap, $r_{\min}$. Its distribution is related to the probability to have a hole of radius $r_{\min}$ at the center, a question which has been studied in [38], but for which more remains to be done in $d > 1$.

Another natural question is how this limiting distribution depends on temperature $T \ge 0$. Previous studies of matrix models related to finite temperature fermions [39, 40] have not addressed questions about extreme value statistics
(see however [36] in connection with 1d fermions). In $d=1$, the distribution of the position of the rightmost fermion in a harmonic trap has recently been studied for finite temperatures [8, 11]. The cumulative distribution, properly centered and scaled, was found to converge to a Gumbel form (i) at very high temperature $T \gg N$ and (ii) also at intermediate temperature where $N^{1/3} \ll T \ll N$. The scale factors associated with these Gumbel forms depend on $T$ and $N$, however these dependencies are rather different in the two regimes. In this work we have investigated what happens in higher dimensions as temperature increases. We have studied only the regime $T \sim N^{1/(3d)}$ and shown that the distribution of $r_{\max}(T)$ is also a Gumbel distribution, although with different parameters to those at $T=0$. The other extreme limit is when $T \to \infty$, where we expect that the fermions behave as independent classical particles in a harmonic potential (the Pauli exclusion principle becomes completely irrelevant). Consequently, the maximal radial distance will correspond to the maximum among a set of $N$ independent and identically distributed random variables, each with a Gaussian tail. Hence we would expect that the limiting cumulative distribution of the maximal distance, properly centered and scaled, will again converge to the standard Gumbel form in all dimensions $d > 1$. How the limiting distribution crosses over from the zero temperature Gumbel form to the infinite temperature Gumbel form as $T$ increases, is an interesting open problem.

Finally, let us stress that some of the predictions of the present theory are in principle testable in cold atom experiments, due to the recent progress in Fermi gas quantum gas microscopes [41–48]. Even though in general the atoms are interacting, the current experimental setups also allow one to tune the strength of the interactions to zero and are able to probe the non-interacting limit [1, 2]. In such microscopes, individual atoms from a correlated many-body system can be imaged in situ with a resolution comparable to the inter-particle spacing, providing direct access to the modulus of the system’s wave-function. A quantum gas microscope gives access not only to the real space distribution but, using a time of flight technique, can also provide the momentum distribution of the dilute gas. Hence both $r_{\max}$, $p_{\max}$ and $x_{\max}$ can be studied. In practice one needs to include the thermal effects, which are not negligible in the experiments. We hope that the present study will stimulate efforts in this direction.

## Appendix A: Higher order terms in the trace expansion

In subsection IIB we used the first term in the trace expansion to derive the distribution of the farthest fermion. Here we consider the behavior of the higher order terms and show that they are indeed negligible in the regime of typical fluctuations. To begin with, we analyze the second term in the trace expansion which is proportional to

$$ \mathrm{Tr}([I_w K_\mu]^2) = \iint_D d\mathbf{x} d\mathbf{y} K_\mu(\mathbf{x}, \mathbf{y})^2. \qquad (A1) $$

where the integration domain $D$ corresponds to $|\mathbf{x}|, |\mathbf{y}| > w$. In general we can write

$$ \int_D dy K_\mu(\mathbf{x}, \mathbf{y})^2 = \int_D dy K_\mu(r\mathbf{n}, \mathbf{y})^2, \qquad (A2) $$

where $r = |\mathbf{x}|$ and $\mathbf{n} = \mathbf{x}/r$. Clearly isotropy implies that Eq. (A2) must be independent of $\mathbf{n}$. We can thus write

$$ \mathrm{Tr}([I_w K]^2) = S_d \int_w^\infty dr r^{d-1} \int_D dy K_\mu(r\mathbf{n}, \mathbf{y})^2, \qquad (A3) $$

where $\mathbf{n}$ is an arbitrary direction. We now write $r = r_{\mathrm{edge}} + a'_n$ and, anticipating that the integral over $\mathbf{y}$ is dominated by points close to $\mathbf{x}$, we also write $\mathbf{y} = r_{\mathrm{edge}}\mathbf{n} + \mathbf{b}' = (r_{\mathrm{edge}} + b'_n)\mathbf{n} + \mathbf{b}'_t$ where $\mathbf{b}'_t$ is the transverse component of the vector $\mathbf{y}$ in the plane whose normal is $\mathbf{n}$. Using the fact that $r_{\mathrm{edge}}$ is large means that we can write

$$ \mathrm{Tr}([I_w K]^2) \approx S_d r_{\mathrm{edge}}^{d-1} \int_{w-r_{\mathrm{edge}}}^{\infty} da'_n \int_{w-r_{\mathrm{edge}}}^{\infty} db'_n \int d\mathbf{b}'_t K([r_{\mathrm{edge}} + a'_n]\mathbf{n}, [r_{\mathrm{edge}} + b'_n]\mathbf{n} + \mathbf{b}'_t)^2. \qquad (A4) $$

In Ref. [8, 11] it was shown that for $\mathbf{x}, \mathbf{y}$ near a point $r_{\mathrm{edge}}\mathbf{n}$, where $\mathbf{n}$ is a unit vector (on the unit sphere), the kernel takes the scaling form

$$ K_\mu(\mathbf{x}, \mathbf{y}) \approx \frac{1}{w_N^d} K_{d, \mathbf{n}}^{\mathrm{edge}} \left( \frac{\mathbf{x} - r_{\mathrm{edge}}\mathbf{n}}{w_N}, \frac{\mathbf{y} - r_{\mathrm{edge}}\mathbf{n}}{w_N} \right), \qquad (A5) $$

where $K_{d, \mathbf{n}}^{\mathrm{edge}}(\mathbf{a}, \mathbf{b})$ is given in (144) in terms of the Airy kernel $K_{\mathrm{Ai}}(a, b) = \int_0^{+\infty} du \mathrm{Ai}(a+u)\mathrm{Ai}(b+u)$ and the integral over the transverse Fourier variable $\mathbf{q}_t$ is over the $(d-1)$-dimensional space. Using the scaling form in Eq. (A5) into
Eq. (A4) and writing $a'_n = w_N a_n$ and $b' = w_N b$ now yields

$$ \mathrm{Tr}([I_w K]^2) \approx S_d \frac{r_{\mathrm{edge}}^{d-1}}{w_N^{d-1}} \int_{\tilde{w}}^{\infty} da_n \int_{\tilde{w}}^{\infty} db_n \int db_t K_{d,n}^{\mathrm{edge}}(a_n \mathbf{n}, b_n \mathbf{n} + \mathbf{b}_t)^2, \quad (A6) $$

where $\tilde{w}$ is as defined in Eq. (38).

We now use the representation Eq. (144) twice to represent $K_{d,n}^{\mathrm{edge}}^2$, with two transverse Fourier variables $\mathbf{q}$ and $\mathbf{q}'$ for each factor. The integral over $b_t$ can be carried out yielding a term $(2\pi)^{d-1}\delta(\mathbf{q}_t + \mathbf{q}'_t)$ and leads to (see also formula (145))

$$ \mathrm{Tr}([I_w K_{\mu}]^2) = S_d \frac{r_{\mathrm{edge}}^{d-1}}{w_N^{d-1}(2\pi)^{d-1}} \int_{\tilde{w}}^{\infty} da_n \int_{\tilde{w}}^{\infty} db_n \int dq_t K_{\mathrm{Ai}}^2(a_n + q_t^2, b_n + q_t^2). \quad (A7) $$

The remaining integral over $\mathbf{q}_t$ can now be carried out and the result written as

$$ \mathrm{Tr}([I_w K_{\mu}]^2) = \frac{S_d S_{d-1}}{2(2\pi)^{d-1}} \left(\frac{r_{\mathrm{edge}}}{w_N}\right)^{d-1} \int_0^{\infty} du \, u^{\frac{d-3}{2}} \int_{\tilde{w}+u}^{\infty} da_n \int_{\tilde{w}+u}^{\infty} db_n K_{\mathrm{Ai}}^2(a_n, b_n). \quad (A8) $$

Using the notation of functional traces we see the above can be simply written as

$$ \mathrm{Tr}([I_w K_{\mu}]^2) = \frac{1}{\Gamma(d-1)} \left(\frac{r_{\mathrm{edge}}}{w_N}\right)^{d-1} \int_0^{\infty} du \, u^{\frac{d-3}{2}} \mathrm{Tr}([P_{\tilde{w}+u} K_{\mathrm{Ai}}]^2), \quad (A9) $$

where $P_s$ is the projector on $[s, +\infty[$. The computation of higher order traces can be done in an identical manner and it is easy to see that

$$ \mathrm{Tr}([I_w K_{\mu}]^n) = \frac{1}{\Gamma(d-1)} \left(\frac{r_{\mathrm{edge}}}{w_N}\right)^{d-1} \int_0^{\infty} du \, u^{\frac{d-3}{2}} \mathrm{Tr}([P_{\tilde{w}+u} K_{\mathrm{Ai}}]^n), \quad (A10) $$

for $n \ge 2$.

In order to analyze the large $\tilde{w}$ behavior, we use the “Christoffel-Darboux formula” for the Airy kernel, in terms of the Airy function $\mathrm{Ai}$ and its derivative $\mathrm{Ai}'$,

$$ K_{\mathrm{Ai}}(x, y) = \frac{\mathrm{Ai}(x)\mathrm{Ai}'(y) - \mathrm{Ai}(y)\mathrm{Ai}'(x)}{x - y}, \quad (A11) $$

and for large $x$ we use the asymptotic formula

$$ \mathrm{Ai}(x) \sim \frac{\exp(-\frac{2}{3}x^{\frac{3}{2}})}{2\sqrt{\pi}x^{\frac{1}{4}}}, \quad (A12) $$

which means that

$$ \mathrm{Ai}'(x) \sim -\sqrt{x}\mathrm{Ai}(x). \quad (A13) $$

Thus, for both $x$ and $y$ large, we have

$$ K_{\mathrm{Ai}}(x, y) \sim \frac{\mathrm{Ai}(x)\mathrm{Ai}(y)}{\sqrt{x} + \sqrt{y}}. \quad (A14) $$

Now consider the second order contribution to the trace formula, specifically the last term of the integrand in Eq. (A9) which can be written as

$$ \mathrm{Tr}([P_{\tilde{w}+u} K_{\mathrm{Ai}}]^2) = \int_0^{\infty} da \int_0^{\infty} db K_{\mathrm{Ai}}^2(a + \tilde{w} + u, b + \tilde{w} + u). \quad (A15) $$

In the above, $\tilde{w}$ is assumed to be large and $u$ is positive so we can use the asymptotic result from Eq. (A14) to obtain

$$ \mathrm{Tr}([P_{\tilde{w}+u} K_{\mathrm{Ai}}]^2) \approx \int_0^{\infty} da \int_0^{\infty} db \frac{\mathrm{Ai}^2(a + \tilde{w} + u)\mathrm{Ai}^2(b + \tilde{w} + u)}{(\sqrt{a + \tilde{w} + u} + \sqrt{b + \tilde{w} + u})^2}. \quad (A16) $$
The asymptotic expansion of the integrand above for large $\tilde{w}$ is dominated by the exponential term in the Airy function and we obtain

$$ \mathrm{Tr} ([P_{\tilde{w}+u} K_{\mathrm{Ai}}]^2) \approx \frac{\exp(-\frac{8}{3}\tilde{w}^{\frac{3}{2}} - 4\tilde{w}^{\frac{1}{2}}u)}{64\pi\tilde{w}^{\frac{5}{2}}}, \qquad (\mathrm{A}17) $$

which yields

$$ \mathrm{Tr}([I_w K_\mu]^2) \sim \left(\frac{r_{\mathrm{edge}}}{w_N}\right)^{d-1} \frac{\exp(-\frac{8}{3}\tilde{w}^2)}{\tilde{w}^{\frac{9+d}{4}}}. \qquad (\mathrm{A}18) $$

In the above if we substitute $\tilde{w} \approx a_\mu = (\frac{d-1}{2} \ln \mu)^{\frac{2}{3}}$ we find that, up to logarithmic corrections

$$ \mathrm{Tr}([I_w K_\mu]^2) \sim \mu^{-\frac{2}{3}(d-1)}. \qquad (\mathrm{A}19) $$

It can be shown by a simple extension of the above argument that

$$ \mathrm{Tr}([I_w K_\mu]^n) \sim \mu^{-\frac{2(n-1)}{3}(d-1)} \qquad (\mathrm{A}20) $$

in the regime of typical fluctuations. Hence, for $d > 1$ the higher order traces are irrelevant in that regime.

[1] I. Bloch, J. Dalibard, W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).

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

[3] G. D. Mahan, Many particle physics, Plenum, NY (1981).

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

[5] W. Kohn, A. E. Mattsson, Phys. Rev. Lett. 81 3487 (1998).

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

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

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

[9] D. S. Dean, P. Le Doussal, S. N. Majumdar, G. Schehr, Europhys. Lett. 112, 60001 (2015).

[10] R. Marino, S. N. Majumdar, G. Schehr, P. Vivo, Phys. Rev. E 94, 032115 (2016).

[11] D. S. Dean, P. Le Doussal, S. N. Majumdar and G. Schehr, preprint arXiv:1609.04366, to appear in Phys. Rev. A (2016).

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

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

[14] Y. V. Fyodorov, Introduction to the Random Matrix Theory: Gaussian Unitary Ensemble and Beyond, Lond. Math. Soc. Lect. Note Ser. 322, 31 (2005).

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

[16] E. J. Gumbel, Statistics of Extremes (Columbia University Press, New York, 1958).

[17] J. Galambos, The Asymptotic Theory of Extreme Order Statistics (R.E. Krieger Publishing Co., Malabar, 1987).

[18] J. Ginibre, J. Math. Phys. 6, 440 (1965).

[19] B. A. Khoruzhenko, H. J. Sommers, Non-Hermitian ensembles, in The Oxford Handbook of Random Matrix Theory, G. Akemann, J. Baik, P. Di Francesco (Eds.), Oxford University Press, Oxford (2011).

[20] B. Rider, D. Sinclair, Ann. Appl. Probab. 4 1621 (2014).

[21] B. Rider, J. Phys. A: Math. Gen. 36, 3401 (2003).

[22] F. D. Cunden, F. Mezzadri, P. Vivo, J. Stat. Phys. 164, 1062 (2016).

[23] D. Chafaï, S. Péché, J. Stat. Phys. 156, 368 (2014).

[24] K. Johansson, Random matrices and determinantal processes, in Lecture Notes of the Les Houches Summer School 2005 (A. Bovier, F. Dunlop, A. van Enter, F. den Hollander, and J. Dalibard, eds.), Elsevier Science, (2006); arXiv:math-ph/0510038.

[25] A. Borodin, Determinantal point processes, in The Oxford Handbook of Random Matrix Theory, G. Akemann, J. Baik, P. Di Francesco (Eds.), Oxford University Press, Oxford (2011).

[26] M. Moshinsky, Y.F. Smirnov, The Harmonic oscillator in modern physics, (Contemporary Concepts in Physics Volume 9) Harwood Academic Publishers Amsterdam (1996).

[27] C. Nadal, S. N. Majumdar, Phys. Rev. E 79, 061117 (2009).

[28] J. Wishart, Biometrica 20, 32 (1928).

[29] P. Vivo, S. N. Majumdar, O. Bohigas, J. Phys. A: Math. Theor. 40, 4317 (2007).

[30] K. Johannsson, Comm. Math. Phys. 209 (2000).
[31] I. M. Johnstone, Ann. Stat. **29**, 295 (2001).

[32] P. J. Forrester, Forum Math. **18**, 711 (2006)

[33] A. Edelman, M. La Croix, Random Matrices: Theory Appl. **04**, 1550021 (2015).

[34] F. Bornemann, M. La Croix, preprint arXiv:1502.05946 (2015).

[35] D. S. Dean, S. N. Majumdar, Phys. Rev. E **77**, 041108 (2008).

[36] K. Johansson, Probab. Theory Rel. **138**, 75 (2007).

[37] P. Le Doussal, S. N. Majumdar, G. Schehr, EPL **113**, 60004 (2016). We use formula (55) and (57) of the Supp. Mat. of the arXiv version (arXiv:1505.01543), with the correspondence $t \to b^3$ and $y \to \tilde{s}$.

[38] S. Torquato, A. Scardicchio, C. E. Zachary, J. Stat. Mech. P11019, (2008).

[39] A. M. Garcia-Garcia, S. M. Nishigaki, J. J. M. Verbaarschot, Phys. Rev. E **66**, 016132 (2001).

[40] A. M. Garcia-Garcia, J. J. M. Verbaarschot, Nucl. Phys. B **586**, 668 (2000).

[41] L. W. Cheuk, M. A. Nichols, M. Okan, T. Gersdorf, R. Vinay, W. Bakr, T. Lompe, M. Zwierlein, Phys. Rev. Lett. **114**, 193001, (2015).

[42] E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaudecerf, G. D. Bruce, S. Kuhr, Nature Physics **11**, 738 (2015).

[43] M. F. Parsons, F. Huber, A. Mazurenko, C. S. Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt, M. Greiner, Phys. Rev. Lett. **114**, 213002 (2015).

[44] W. S. Bakr, J. I. Gillen, A. Peng, S. Folling, M. Greiner, Nature **462**, 74 (2009).

[45] J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, S. Kuhr, Nature **467**, 68 (2010).

[46] W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Fölling, L. Pollet, M. Greiner, Science **329**, 547 (2010).

[47] A. Omran, M. Boll, T. A. Hilker, K. Kleinlein, G. Salomon, I. Bloch, Ch. Gross, Phys. Rev. Lett. **115**, 263001 (2015).

[48] G. J. A. Edge, R. Anderson, D. Jervis, D. C. McKay, R. Day, S. Trotzky, J. H. Thywissen, Phys. Rev. A **92** 063406 (2015).