Full counting statistics for interacting trapped fermions

We study $N$ spinless fermions in their ground state confined by an external potential in one dimension with long range interactions of the general Calogero-Sutherland type. For some choices of the potential this system maps to standard random matrix ensembles for general values of the Dyson index $\beta$. In the fermion model $\beta$ controls the strength of the interaction, $\beta=2$ corresponding to the noninteracting case. We study the quantum fluctuations of the number of fermions ${\cal N}_{\cal D}$ in a domain $\cal{D}$ of macroscopic size in the bulk of the Fermi gas. We predict that for general $\beta$ the variance of ${\cal N}_{\cal D}$ grows as $A_{\beta} \log N + B_{\beta}$ for $N \gg 1$ and we obtain a formula for $A_\beta$ and $B_\beta$. This is based on an explicit calculation for $\beta\in\left\{ 1,2,4\right\} $ and on a conjecture that we formulate for general $\beta$. This conjecture further allows us to obtain a universal formula for the higher cumulants of ${\cal N}_{\cal D}$. Our results for the variance in the microscopic regime are found to be consistent with the predictions of the Luttinger liquid theory with parameter $K = 2/\beta$, and allow to go beyond. In addition we present families of interacting fermion models in one dimension which, in their ground states, can be mapped onto random matrix models. We obtain the mean fermion density for these models for general interaction parameter $\beta$. In some cases the fermion density exhibits interesting transitions, for example we obtain a noninteracting fermion formulation of the Gross-Witten-Wadia model.


Overview
The full counting statistics (FCS), which measures the fluctuations of the number of particles N D inside a domain D has been studied extensively in the context of shot noise [1], quantum transport [2], quantum dots [4,5], non-equilibrium Luttinger liquids [3] as well as in quantum spin chains and fermionic chains [6][7][8][9][10][11]. The FCS is particularly important for noninteracting fermions because of its connection to the entanglement entropy [12][13][14][15][16]. For free fermions, in the absence of external potential and at zero temperature it is well known that both the variance of N D and the entropy grow as ∼ R d−1 log R with the typical size R of the domain D in space dimension d [18][19][20][21][22][23].
Recently these results have been extended for noninteracting fermions in the presence of a confining potential. This is important for applications e.g. to cold atoms experiments [26][27][28][29] where the fermions are in traps of tunable shapes [30,31]. In a confining potential, the Fermi gas is supported over a finite domain. Its mean density is inhomogeneous and can be calculated using the well known local density approximation (LDA) [24,25]. To compute quantum correlations, in particular at the edge of the Fermi gas, where the density vanishes and the LDA method fails, more elaborate methods have been developed [32][33][34].
In d = 1, one can exploit the fact that for specific potentials, the problem at zero temperature can be mapped to standard random matrix ensembles, for which powerful mathematical tools are available. For instance, for N noninteracting spinless fermions in a harmonic well, described by the single particle Hamiltonian H = p 2 2 + V (x) with V (x) = 1 2 x 2 , the quantum joint probability distribution function (PDF) for the positions x = {x i } i=1,. ..,N of the fermions takes the form with β = 2 and Ψ 0 ( x) denotes the ground state many body wave function. Remarkably, Eq.
(1), specialised to β = 2, coincides with the joint PDF of the eigenvalues λ i of a random matrix belonging to the Gaussian unitary ensemble (GUE). The two problems thus map into each other with the identification λ i = x i . As a result, for large N , the mean fermion density, i.e., the quantum average ρ(x) = i δ(x − x i ) , takes the Wigner semi-circle form ρ(x) ρ bulk (x) = 1 π √ 2N − x 2 in the bulk, i.e., for x ∈ [x − e , x + e ] and vanishes beyond the edges x ± e ± √ 2N . To discuss the FCS within an interval [a, b] at large N , i.e., the fluctuations of N [a,b] , one needs to distinguish two natural length scales: the microscopic scale given by the interparticle distance ∼ 1/ √ N , and the macroscopic scale of order x + e − x − e ∼ √ N . It is well known since Dyson and Mehta [35,36] that for an interval of microscopic size the variance for the GUE is given by [6,[37][38][39][40][41][42][43][44] Var N [a,b] 1 π 2 log 2N −a 2 |b − a| + c 2 , for √ N |b − a| = O(1) 1, with c 2 = γ E + 1 + log 2, where γ E is Euler's constant. This result is obtained from the celebrated sine kernel which describes the eigenvalue correlations in the GUE at microscopic scales. The formula (2) thus carries to the fermions in the harmonic potential. The FCS for an interval of macroscopic size has been studied more recently. For the harmonic well, some results for the variance in that regime were obtained using the connection to the GUE, both in mathematics [45][46][47][48], and in physics using a Coulomb gas method [41,42].
The mapping to random matrix theory (RMT) holds however only for a few specific potentials. For instance the so-called Wishart-Laguerre ensemble is related to the potential x 2 for x > 0. For an arbitrary smooth V (x), not necessarily related to RMT, we have recently obtained a general formula [49] for the variance Var N [a,b] for a macroscopic interval [a, b] in the bulk. The method used combines determinantal point processes with semi-classical (WKB) approaches. In the special cases related to RMT, the formula recovers the available exact results [50]. It is also in general agreement with recent approaches relying on inhomogeneous bosonization [51][52][53], although our method allows for more precise and controled results.
One can also ask about higher cumulants of N [a,b] , i.e., beyond the variance given in (2), both for microscopic and macroscopic interval [a, b]. In the absence of potential, i.e., for free fermions, there exist results for the higher cumulants which are obtained using the sine-kernel [6,40,44]. A natural conjecture, which we put forward in [49] is that these higher cumulants are determined solely from fluctuations on microscopic scales. Consequently (i) they are independent of the size of the interval (within the bulk) and (ii) they are universal, i.e., independent of the precise shape of the potential (assumed to be smooth). This conjecture was used to obtain a prediction for the entanglement entropy for noninteracting fermions in a potential in [49].
An outstanding question is the study of the FCS for interacting particles [54]. It is of current interest for cold atom experiments, which have recently measured particle number fluctuations in the 1d Bose gas [55,56]. On the theory side however there are only a few results, even for integrable models. For instance, for the delta Bose gas (Lieb-Liniger) model, an exact formula was derived using the Bethe ansatz for the FCS in the limit of a very small interval [57,58]. Results for larger intervals were also obtained, but are only valid in the weak interaction/high temperature regime [59]. For interacting fermions, numerical results were obtained for the Hubbard model [60]. In the context of spin chains, several exact formulae for the FCS were obtained, e.g for the XXZ spin chain [61,62], for the transverse field Ising model [10] and for the Haldane-Shastry chain [9].
In view of our previous works on the FCS of noninteracting trapped fermions it is thus natural to look for extensions which include interactions. A promising direction is to explore further the connection between random matrix theory for a general value of the Dyson index β, and trapped fermions in 1d in the presence of two-body interactions of the Calogero-Sutherland type [63][64][65][66]68]. The simplest example corresponds to the Gaussian β ensemble (GβE) which contains the GUE for β = 2, as well as the other standard ensembles, the GOE for β = 1 and the GSE for β = 4 [39,68]. For general β they can be constructed from random tridiagonal matrices [69]. The joint PDF of their eigenvalues λ i is given by (1) with the substitution x i ≡ β 2 λ i . The important observation is that Eq. (1) is also the quantum joint PDF, |Ψ 0 ( x)| 2 , of the positions x i of N fermions in the ground state of the following N body Hamiltonian It thus describes fermions which interact through either repulsive (β > 2) or attractive (1 ≤ β < 2) long range 1/x 2 interaction. As we discuss below there are several other examples of two-body interactions and trapping potentials for which a similar connection exists. Note that there are also lattice versions of these models, e.g. Haldane-Shastry spin chains corresponding to a discretized version of the circular β ensemble (CβE) with β = 4, for which FCS results exist [9]. In this paper we will use the relations between interacting fermions and RMT for general β, to obtain precise predictions for the FCS for various examples of trapped fermions in the presence of interactions. In the rest of this section we first present the main models that we will study, we explain the main idea of the method and we present the main results.

Fermions' domain
Fermion potential V (x) Fermion interaction W (x, y)

RMT ensemble
Matrix potential V 0 (λ) Map λ(x) 16 sin 2 π(x−y)  Table 1: The mappings between (i) models of interacting trapped fermions studied here and (ii) the standard random matrix ensembles. The variable x denotes the positions of the fermions, λ the eigenvalues of the RMT ensemble, and the mapping λ(x) is displayed in the last column. The first three columns denote respectively the domain for the fermions, the external potential V (x), and their interaction W (x, y), defined in (3). Note that in the second line periodic boundary conditions are to be understood for the fermionic system. The next two columns indicate the RMT ensemble and the matrix potential V 0 (λ), see Eq. (6). Here β is the Dyson index which varies continuously and corresponds to noninteracting fermions for β = 2.

Models and mappings
Let us now describe the class of models which we study in this paper. In this section we focus on models related to RMT, while further extensions will be discussed below. Here we consider N spinless fermions trapped in an external potential V (x) and with two-body interactions parameterized by a symmetric function W (x, y) = W (y, x). The Hamiltonian is of the general form (we use units such that m = = 1) In this paper we study specific choices for V (x) and W (x, y) such that the joint PDF of the positions of the fermions in the ground state can be written in the form with w a symmetric function w(x, y) = w(y, x). One example of such models corresponds to Eqs.
Let us briefly review the ensembles of interest. We denote λ i , i = 1, . . . , N the eigenvalues of a random matrix in an ensemble such that the joint PDF can be written as ...,N and Z N is a normalisation constant. Here β is the Dyson index and V 0 the matrix potential (not to be confused with the fermion potential V ). For instance the Gaussian-beta ensemble (GβE) corresponds to the case where the eigenvalues are on the real axis, λ i ∈ R, with V 0 (λ) = β 2 λ 2 . It contains the Gaussian unitary, orthogonal and symplectic ensemble for β = 2, 1, 4 respectively. The circular-beta ensemble (CβE) corresponds to the case where the eigenvalues are on the unit circle in the complex plane, with V 0 (λ) = 0, and includes the circular unitary ensemble (CUE) for β = 2. The Wishart-Laguerre-beta ensemble (WLβE) corresponds to λ i ∈ R + and V 0 (λ) = β 2 λ − γ log λ. The Jacobi-beta ensemble (JβE) corresponds to λ i ∈ [0, 1] and V 0 (λ) = −γ 1 log λ − γ 2 log(1 − λ). In all these cases and for any β, Z N has an explicit expression as a Selberg integral [68], and the ensembles can be mapped onto certain tridiagonal matrices [69]. These ensembles are recapitulated in the Table 1.
The general idea behind the mapping between fermions and RMT is that, upon some map which we denote λ(x), i.e., λ i = λ(x i ), one can identify the joint PDF (6) with the quantum joint PDF (5) corresponding to the many-body ground state of the fermion system with Hamiltonian (4). Taking into account the Jacobian of the map, the correspondence reads The simplest case is the mapping λ(x) = 2 β x from the GβE to the fermions on the real axis described by the model (3). For the CβE ensemble the map is λ(x) = e ix 2π L , where the fermions live on the periodic ring, x j ∈ [0, L]. In that case it maps onto the Sutherland model, without any external potential V (x) = 0, see second line of the Table 1. For the WLβE ensemble the map is λ(x) = 2 β x 2 and the fermions live on R + , x i > 0, with potential sum of harmonic and 1/x 2 wall, and 1/x 2 type interactions as given in the third line of the Table 1. Finally for the JβE ensemble the map is λ(x) = 1 2 (1 − cos x) and the fermions live in a box x i ∈ [0, π] with potential and interactions given in the fourth line of the Table 1. For γ 1 = γ 2 = 1/2 the potential is a hard box with Dirichlet boundary conditions. For β = 2 the fermions are noninteracting in all four cases, and their positions x i in the ground state form a determinantal point process. Note that the above mappings are valid for any N .
To summarize, our main strategy behind the mapping between the ground state of trapped fermions with two-body interactions W and the joint distribution of the eigenvalues of a matrix model consists of the following three steps.
• We consider the Hamiltonian H N in (4) which has only one and two-body potentials, V and W respectively. We then write the many-body ground state wave function in any given ordered sector, e.g.
is of the form (5) consisting only of one-body and two-body terms.
• We next substitute this wave function Ψ 0 ( x) ∼ e −U ( x)/2 in the Schrödinger equation H N Ψ 0 = E 0 Ψ 0 (in an ordered sector). The main condition is that this equation is satisfied for some value of the ground state energy E 0 , i.e., that no three-body interaction is generated upon applying the kinetic operator. This condition selects some special families of potentials V and interactions W . In the absence of potential this approach dates back to Sutherland and Calogero [66,67]. This is a standard although tedious calculation, recalled in Appendix A, allowing also to determine E 0 for each model. The fact that Ψ 0 is indeed the ground state is ensured by the additional condition that Ψ 0 ( x) vanishes only at x i = x j for i = j, but not elsewhere [64], see also [63,71].
• Finally, we identify the quantum probability, given by |Ψ 0 ( x)| 2 , as the joint PDF of the eigenvalues λ 1 , . . . , λ N of a random matrix, under a map λ i = λ(x i ), and we show how to construct this map explicitly for several examples. This last step allows us to identify new connections between interacting (and noninteracting) fermions and random matrix models, see e.g. Section 6.
Mean density. The simplest observable to compute is the average density ρ(x) of the fermions where . . . denotes expectation values with respect to the ground state. In particular one can ask how the interactions modify this density as compared to the noninteracting case W = 0. In the large N limit and in the absence of interactions, the density in the bulk reads ρ(x) 1 π 2(µ − V (x)) + as given by LDA or semi-classical methods (we denote everywhere (x) + = max(x, 0)). Note that the LDA works only for noninteracting fermions and in the bulk [34]. Here µ denotes the Fermi energy which is determined by the normalization condition dxρ(x) = N (e.g., µ N for the harmonic oscillator (HO) considered above). Note that for some integrable systems, the LDA may be improved as in Ref. [108] to include interactions. We will not explore this route here.
For the models in Table 1, the noninteracting case corresponds to β = 2. To obtain the density for arbitrary β one can interpret the PDF (5), or equivalently (6), as the Boltzmann distribution for a gas of classical particles at unit temperature, with energy U ( x), or equivalently F ( λ). Using (7) and (8), in both cases, the interaction between these particles is logarithmic which corresponds to the 2d Coulomb interaction. In the large N limit and in the presence of a confining potential, the equilibrium density is obtained by minimizing the corresponding energy. This Coulomb gas (CG) method has been widely used in the context of RMT. Rewriting (6) one immediately sees that the Coulomb gas result for a general β coincides with that of a gas with β = 2 and a matrix potential 2V 0 (λ) /β. Using the known results for the average eigenvalue density, defined as σ(λ) = 1 N i δ(λ − λ i ) , we can write, respectively for the GβE and WLβE The subscripts 'W' and 'MP' stand for Wigner (semi-circle) and Marcenko-Pastur, respectively. Using the mapping to the fermions with and λ(x) = 2 β x and λ(x) = 2 β x 2 for the GβE and WLβE respectively, we obtain the fermion density for the models in the first and third line of the Table 1 as where θ(x) is the Heavisde function. The positions of the two edges are thus x = x ± e ± √ βN in the first case, while x − e 0 and x + e √ 2βN in the second one. For β = 2 it agrees with the LDA result, and it shows that the Fermi gas expands for β > 2 and shrinks for β < 2, as compared to the noninteracting case β = 2, while retaining a semi-circular shape. In the case of the box, corresponding to the JβE, the density is uniform in the large N limit. Note that in the large N limit, with γ = O(1), the 1/x 2 part of the potential in the third and fourth models in the Table 1 do not affect the bulk density (for a different scaling see below). They become important only in the region close to the wall (see below).

Outline and main results
In this paper we study the statistics of N (β) I , i.e., the number of fermions in an interval I, for the models in the Table 1 for any β ≥ 1, and, in a second stage, for a larger class of models. In parallel to the applications to fermions we also obtain new results in the corresponding random matrix ensembles, with a slightly larger domain of validity, i.e., for any β > 0.
In Section 2 we study the variance of N (β) [a,b] for an interval I = [a, b] of macroscopic size in the bulk for β = 1, 2, 4 in the large N limit, and then propose an extension to any β. In all these cases the variance grows logarithmically with N for large N , and we obtain the amplitude of the logarithm together with the O(1) correction term which has a non trivial dependence on the two edges a, b on macroscopic scales. In the noninteracting case β = 2 there exists a formula, recalled here in Eq. (24) for the variance VarN for a general potential V (x). For the harmonic potential V (x) = x 2 2 , which corresponds to GβE, we extend this formula to β = 1, 2, 4, and it reads βN are the positions of the two edges, as can be seen in Eq. (13). For β = 1, 2, 4 the constant c β takes the values Here we argue that formula (14) extends to the model (3) of interacting fermions with general β in the harmonic potential. Using related works [72,73] (see discussion below) we propose the following expression as a series representation for c β where here and below ψ (k) (z) = d k+1 dz k+1 log Γ(z) is the polygamma function. We have checked numerically the predictions (14), (17) for the variance (see Fig. 4).
In fact, going beyond the harmonic potential, our more general prediction for the models in Table 1 reads at large N and in the bulk [74] βπ 2 2 VarN where a = a 2/β and b = b 2/β for the models in line 1 and 3 in the Table 1 and a = a and b = b for the other two models (on a circle and in a box). On the right hand side of Eq.
is the variance for noninteracting fermions (i.e., for β = 2) in the presence of a potential V (x) indicated in the Table 1 and given by the general formula (24) (which takes simpler forms for the models in the Table 1). The constant c β is independent of the model, and (18) also holds on microscopic scales in the limit of large interval (as in Eq. (2)). Finally, in Section 2 we also analyze the case of a "semi-infinite" interval, i.e., [a, +∞[ for the GβE, [0, b] for WLβE and JβE, for which we have a similar prediction. One consequence of our prediction (18) is that in the microscopic limit (|a − b| small compared to the size of the Fermi gas) one has [75] Var for k F (a)|b − a| = O(1) 1. In Section 3 we study the higher cumulants of N [a,b] . We present the following conjecture for interacting fermions. Consider an interval [a, b] inside the bulk. For the models displayed in Table 1 and for any β, the cumulants of N (β) [a,b] of order 3 and higher are determined solely from the microscopic scales. This implies that these higher cumulants are identical to those of the CβE. The cumulants for the CβE have been given in [72], using yet another conjecture about extended Fisher-Hartwig asymptotics for CβE, formulated in [73]. We will thus use these formulae and obtain here the full counting statistics for a larger class of interacting fermion models. These cumulants for general β admit the following series representations for arbitrary integer p > 1, while the odd cumulants vanish. For β ∈ {1, 2, 4} the explicit evaluation for the fourth cumulant gives where ζ(z) is the Riemann-zeta function. The conjecture extends naturally to the case of an interval with only one point in the bulk (i.e., for a "semi-infinite" interval). It is a natural extension of the conjecture previously formulated in Ref. [49] for noninteracting fermions β = 2 and recalled in the introduction. One can check that formula (20) reduces to Eq. (23) in [49] in the case β = 2. This conjecture can be checked in a few cases, with impressive agreement. For instance in Section 4 we study the limit from the bulk to the edge for any β. In the case β ∈ {1, 2, 4} we can compare with the results of Bothner and Buckingham [77] obtained by Riemann-Hilbert methods. We find that it agrees in a quite non-trivial way.
In Section 5 we discuss the approaches to interacting fermions using bosonisation in terms of the Luttinger liquid. As we explain, the Luttinger parameter is given here by K = 2/β for the models in Table 1.
Finally in Section 6 we present a more general class of interacting fermion models which have a ground state wave function of the one-and two-body form (5). Some of these models still map onto random matrices, with however a more general matrix potential V 0 . We study in detail the example of fermions on the circle in the presence of an external periodic potential which, in the noninteracting case β = 2 turns out to be related to the so-called Gross-Witten-Wadia model in high energy physics [78,79]. The density in this model exhibits an interesting transition, and we show that the LDA formula in that case reproduces the well known results obtained in Refs. [78,79] from the Coulomb gas method. As discussed there we expect that our results for the counting statistics extend to this more general class of models.
2 Number variance 2.1 Previous results for noninteracting fermions (β = 2) in an external potential In a recent work [49] we have calculated the variance of the number of fermions in a domain D in d = 1, for noninteracting fermions in their ground state in a general potential V (x). In this case the positions of the fermions form a determinantal point process. This means that the n-point correlation function can be written as a n × n determinant built from the so-called kernel K µ (x, y). As a result the variance can be computed from the following formula [39,68] in terms of the kernel. By plugging the large-N asymptotic form of the kernel (given by the WKB expansion of the eigenstates of the single-particle Hamiltonian) into (23) we obtained the leading-and subleading-order terms in the number variance, which are generically of order O(log N ) and O(1) respectively. Consider a confining potential, such that the bulk density For an interval [a, b] in the bulk with |a − b| 1/k F (a), we obtained that for N 1 (i.e., µ 1) the variance is given by [49] where , and c 2 is given in (15). For a semi-infinite interval, and for any a in the bulk, the variance reads [49] VarN For the harmonic potential V (x) = 1 2 x 2 this gives the explicit expression for the variance for the semi-infinite interval [49] VarN [a,+∞ . For an interval in the bulk, Eqs. (24) and (25) lead to the formula given in (14) with β = 2. Note that one can eliminate the Fermi energy µ in all above formula and express all quantities as a function of the number of fermions N using the relation valid at large N . Since the Fermi energy does not have a direct meaning for interacting fermions, it is indeed more natural to use N , in order to study the dependence in β, at fixed N , as we do below. Note that in the interacting case, the zero-temperature chemical potential can be obtained in the large N limit as µ = ∂ N E 0 where E 0 = E 0 (N, β) is the ground state energy (which, in the models studied here can be calculated exactly, see Appendix A). For β = 2 this definition of µ coincides with the Fermi energy.
We will now compute the number variance for interacting fermions, β = 2. For this purpose we recall the definition of a more general observable, the covariance function.

Two-point covariance function
It is useful to define the two point covariance function C (x, y) which gives the covariance between the numbers of particles in infinitesimal intervals around two distinct points x and y: Cov Using the linearity of the covariance, one immediately obtains, for any two nonintersecting domains D 1 and D 2 Using that N = N D + ND, where N is the total number of fermions, and whereD is the complement of D, we obtain a convenient expression for the number variance in a domain: For noninteracting fermions, the determinantal structure can be used in order to express the covariance function in terms of the kernel, C (x, y) = −K µ (x, y) 2 and one recovers the formula in (23).

Number variance for interacting fermions in a harmonic trap
Let us now discuss the case of interacting fermions described by the model (3), which corresponds to random matrices in the GβE. For interacting fermions, the positions of the fermions do not form a determinantal point process, however for β = 1, 4 they exhibit a Pfaffian structure which allows for (complicated) exact expressions for C(x, y) for any finite N (see e.g. [39,68,80,81]). Here, to calculate the variance for β = 1, 2, 4, we will only need their large N asymptotics. We will first recall their expressions separately for the case of microscopic scales |x − y| = O(1/k F (x)) and macroscopic scales |x − y| 1/k F (x). Indeed, when computing the integral in (31), as we do below, there are contributions from both regimes of scales.
(i) Microscopic scales. For x and y in the bulk with x − y microscopic, C(x, y) can be obtained from the Eqs. (18)- (20) in [82] using the mapping from the Gaussian ensembles to the fermions in the harmonic potential (the same formula also holds, see Eqs. (18) and (19) in [82], for the circular ensembles, i.e., for fermions on the circle) where ρ(x) is the fermion density given in (13), and s (r) = sin (πr) πr , Ds (r) = ds dr = πr cos (πr) − sin (πr) Is (r) = r 0 s r dr = Si (πr) π , Js (r) = Is (r) − sgn (r) 2 , where sgn (r) is the sign function and Si (z) = z 0 sin t t dt is the sine integral. Note that for β = 2 the positions of the fermions form a determinantal process, with an associated kernel given in the bulk by the function s(r) which is the well-known sine-kernel.
(ii) Macroscopic scales. For x and y well separated in the bulk, and for the GβE for arbitrary β the covariance function C(x, y) is also known in the large N limit [83][84][85][86]. Under the RMT to fermion mapping it leads to up to rapidly oscillating terms that average out to zero when integrating over macroscopic domains. For noninteracting fermions β = 2, Eq. (38) was also derived directly from the fermion model, see the Supp. Mat. of [49] (for a numerical check of this formula see [87]). Using a Coulomb gas method Eq. (38) was extended to arbitrary matrix potentials and β in [84]. One can check that the above formula match between microscopic and macroscopic scales, i.e., that the limit r 1 in (32) agrees with the limit |x − y| √ βN in (38). The double integral (31) can then be calculated in the large-N limit. The calculation is performed in the Appendix C. First one can approximate C(x, y) 0 if either x or y are not in the bulk. Plugging in C(x, y) from (32) for x near y, and (38) for x far from y (this procedure works because there is a joint regime where both of the approximate expressions for C(x, y) are valid) one finds the result for the variance given in (14) for β ∈ {1, 2, 4}. For instance, for a finite interval [−a, a] centered around the origin and contained in the bulk, a = a/ √ β N < 1, Eq. (14) simplifies into The leading term agrees with a Coulomb gas calculation for general β, [41,42]. In the limit |b −ã| 1, Eq. (14) matches with the microscopic result (2). We also obtain the result for a semi-infinite interval whose edge is in the bulk, |ã| < 1 as In both formulae (39) and (40) the constants c β for β = 1, 2, 4 are related to the so-called Dyson-Mehta constants and given in (15). For β = 1 and β = 4, we have checked numerically the predictions given in (39) and (40) for the fermion model using the correspondence in the first line of the Table 1. We have performed exact diagonalizations of the GβE with β = 1, 4 in Figs. 1 and 2 respectively (the case β = 2 has been tested numerically in [49]). The convergence at large N appears to be slower as β is increased, which is known to occur quite generally, see e.g. Fig. 3

Conjecture for general β and results for other potentials
In the previous section, we have calculated the number variance for interacting fermions in the harmonic potential for β = 1, 2, 4. We now conjecture that the formula (14), (39) and (40) hold for general values of β, with an a priori unknown β-dependent constant c β . The rationale behind this conjecture is that (i) the expression (38) for C(x, y) on macroscopic scale is valid for arbitrary β, a result which comes naturally from the Coulomb gas calculations [41,42,84] (ii) the above calculations for β = 1, 2, 4 show that the β-dependence of the constant part, c β , is determined from microscopic scales only. Hence we expect that it is independent of the RMT ensemble in the Table 1. As we argue below, this constant is given for general β by formula (17), which we will justify in Section 3 based on previous works on the CβE. The conjecture can be expressed as follows. For an arbitrary interval [a, b] in the bulk and for the fermion models listed in Table 1, there is a relation between the variance for an arbitrary β ≥ 1 and the variance for β = 2 (up to a possible rescaling of lengths). This relation is given in (18) above.
Let us now present a few results which follow from this conjecture. For fermions on the circle with L = 2π this leads to In the microscopic limit, this formula agrees with (19) with k F = πρ = N/2. Consider now interacting fermions related to the WLβE in the potential (line 3 in Table  1) In the introduction, we have discussed this model when the parameter γ = O(1) in which case the 1/x 2 hard wall potential does not affect the bulk properties for x > 0, such as the density given in Eq. (13). Another interesting limit amounts to scale the parameter γ ∼ µ ∼ N in which case the effect of the 1/x 2 potential is to open a gap in the density of fermions near the origin. Let us recall that the associated Wishart-Laguerre (WL) matrix potential is V 0 (λ) = β 2 λ − γ log λ. By a similar argument as in Eq. (10) and below, one can absorb the β dependence in the product 2 β V 0 (λ), by rescaling γ. As a result the eigenvalue density σ(λ; γ) for the WLβE takes the scaling form at large N (as obtained e.g. from the Coulomb gas method)   Table 1 associated to the WLβE with γ = 2, and β = 1 (a) and β = 2 (b). The blue markers are the empirical variance computed over 5 × 10 4 simulated WLβE matrices with N = 100, and the red lines are our theoretical prediction (48).
One can now use our main conjecture (18) and its analog for semi-infinite intervals, and the result that we obtained in [49] for the case β = 2, to predict the variance of the number of fermions in an interval for general β for the potential (42). For the interval [0, a] with a in the bulk, this leads to 2βN (which is the position of the edge atγ = 0). The theoretical prediction (48) is compared with numerical simulations of Wishart matrices in Fig. 3 with γ = 2 and β ∈ {1, 2}, with excellent agreement. A formula analogous to (48) for a general interval [a, b] with a, b in the bulk is given in Appendix D, where we also give some formula for the Jacobi box potential (line 4 in Table 1) as well as the details of the derivation of (48).

Higher cumulants
We now study the higher cumulants (larger than 2) of the number of fermions in an interval for the interacting fermion models displayed in Table 1. In our previous work for noninteracting fermions in [49] we had conjectured, and checked with available rigorous results for several potentials, that the higher cumulants are determined solely from microscopic scale. Hence they are independent of the potential in the large N limit. Here we will go one step further and conjecture that this remains true in the interacting case for general β. Although the numerical values of these cumulants depend non trivially on β, i.e., on the interaction strength, they are insensitive to the details of an external smooth potential. Indeed these cumulants are determined at microscopic scales where the 1/x 2 interactions dominate over the local variations of the potential. Consequently we can conjecture that the higher cumulants are the same as for fermions on a circle without a potential, i.e for the CβE. It turns out that the cumulants for the CβE were recently predicted in Ref. [72] in a different context. Consider the periodic model in the second line of Table 1 where x is the coordinate along a circle of perimeter L = 2π. The result of [72] gives the FCS generating function for N [a,b] , i.e., the number of fermions with positions x i ∈ [a, b] as [89] log e up to terms that vanish in the large N limit. Here t is a parameter [90], and for β = 2s/r, with s, r integers mutually prime where the coefficientsC It is obvious from this formula thatC (β) 2p+1 = 0 hence all the odd cumulants vanish, i.e., N 2p+1 [a,b] c = 0. We thus focus now on the even cumulants. Although the coefficientsC (β) k are defined here for rational values of β, it is possible to obtain expressions for these coefficient for any real β, as an explicitly continuous function of β. This is achieved using the fact that any real β can be reached by a sequence β = 2s n /r n of arbitrary large s n , r n and performing an asymptotic analysis (see details in [72]). The result for k = 2p with p ≥ 2 can be written in terms of the following double series In addition, one of the sums (either over ν or over q) can be carried out, leading to two equivalent "dual" expressions is the polygamma function. The above series are convergent for p ≥ 2, since at large x one has ψ (2p−1) (z) The asymptotics for small and large β can be obtained from either of the dual series in (54) and can be found in [72]. For the classical values β ∈ {1, 2, 4} these series can be performed explicitly, e.g. see formula (22) for the fourth cumulant. Note that the formula (53) transforms simply under the "duality" β → 4/β. This duality was studied in [93].
From our conjecture, the formula (51) for the cumulants of the fermion model on the circle (i.e., the CβE) is thus predicted to hold for all the fermion models in Table 1, with no modification. Indeed the rescaling of lengths is unimportant here since the values of these cumulants are independent of the size of the intervals (assumed here to be macroscopic in the bulk). In the case of a semi-infinite interval, e.g. [a, +∞[ for the quadratic potential, the result is divided by a factor of 2.
We can now return to the question of the O(1) term in the second cumulant (the variance) as discussed in the previous section. In particular the above predictions based on the CβE allow to obtain explicitly the universal constant c β which enters in all the formulae for the variance of the fermion models considered here. Comparing the formula (41) with the O(t 2 ) term in (49) one finds that the relation (18) holds, together with c β = log 2 + 1 is given in (52). Using the analysis ofC in [72], the constant c β can be written in several alternative forms, either as a convergent double series or, performing one of the sums, as a convergent simple series as given in the Introduction, see (17), or as the dual series We have tested the prediction for c β (17) numerically, together with the prediction for the variance (40), see Fig. 4. Using the correspondence in Table (1) we have diagonalized GβE matrices generated using the Dimitriu-Edelman tridiagonal matrices [69] for various sizes N .
A rather large even-odd finite N effect is observed, however the average value over consecutive N 's is very close to the predictions. Note that in Fig. 4 we tested these predictions also in the range 0 < β < 1, which is only relevant for RMT (or log-gases) but not for the fermion models studied in the rest of this paper (since in the latter models β ≥ 1). It would be interesting to test our predictions for the higher cumulants too. This is more computationally demanding, but nevertheless the fourth cumulant was tested numerically in [70] for the GUE.

FCS near the edge and matching with the bulk
Until now we have discussed the counting statistics in the large N limit for an interval which has at least one point inside the bulk. Let us consider a general smooth potential V (x) such that the Fermi gas has two edges x x ± e , where the LDA density vanishes. There is a region near these edges, of width denoted w N , where it is known that the quantum fluctuations are enhanced, and that the counting statistics are different from the bulk. While this region has been studied in the noninteracting case (β = 2), there are only a few recent results for the counting statistics in the edge region for the interacting case. They were obtained in the context of RMT, specifically for the GβE for β = 1, 4. This corresponds to interacting fermions in a quadratic potential. In Ref. [77] the FCS (i.e., all the cumulants) have been obtained, in the outer edge region, i.e., for an interval [a, +∞) where x + −a w N is O(1) but large (i.e., in the crossover region from the edge to the bulk). Inside the edge region, there are recent results about the second cumulant for general linear statistics for β = 1, 4 [95]. In this section we discuss how these results compare with our conjecture for the cumulants for general β. We start by recalling the noninteracting case and the matching between the bulk and the edge regions.

Noninteracting case β = 2
In the case of noninteracting fermions the positions x i form a determinantal point process based on the kernel K µ (x, y) discussed in Section 2. For any smooth confining potential V (x), as discussed in [34], for x, y near the right edge x + (and similarly for x − ) the kernel takes the universal scaling form K µ (x, y) Here K Ai is the Airy kernel given by  (17). One observes a (rather) slow convergence of the numerical results to our conjecture as N is increased. In (c) and (d), analogous results are plotted for 0 < β < 1 (and N = 2000, 2001), again with good agreement between the numerical results and our conjecture. These results are only relevant for RMT (or log-gases) but not for the fermion models studied in the rest of this paper (since in the latter models β ≥ 1).
Using this scaling form one obtains the number variance for any interval in the edge region in terms of the Airy kernel as was done in [41,42] for the case of the harmonic oscillator/GUE, where the scaling function V 2 (â), defined in [41,42], is universal, i.e., independent of the potential V (x), in terms of the scaling variableâ. One interesting question is the matching of the number variance as a in (58) moves from the bulk to the edge. In [49] it was shown that the asymptotic behavior for a → x + coming from the bulk reads in terms of the edge scaling variableâ defined in (58). The result (59) matches exactly with the formula (58) in the limitâ → −∞, which corresponds to the crossover from the edge to the bulk. The comparison between the two results is performed in Appendix F. This crossover was also obtained in the case of the harmonic potential (i.e., for the GUE) in [77] (see also [96]). In fact in Ref. [77] the FCS, i.e., the higher cumulants were also given for β = 2 in this crossover regime. As we discuss below and in Appendix E, this result matches perfectly our predictions for the higher cumulants in the bulk given in [49].

Interacting case
Let us start with the case of the harmonic potential V (x) = 1 2 x 2 . From the RMT-fermion correspondence in the first line of Table 1, i.e., λ(x) = 2/β x, one finds that for general β the right edge is at position x + = √ βN and the width of the edge region is It is known in RMT that the edge properties of the GβE are described by the Airy β point process denoted a β i , which implies that the fermion positions in the edge region and for N → +∞ can be written as The statistics of the a β i is described by the so-called stochastic Airy operator [97][98][99]. It is known that the largest eigenvalue (i.e., the rightmost fermion) is described by the β-Tracy-Widom distribution Prob(max i a β i < s) = F β (s) which depends continuously on β. Explicit expressions in terms of Fredholm determinants or solutions to Painlevé equations are known for β ∈ {1, 2, 4} [100,101] (see also [77,102,103]). From known results for the mean density σ(λ) of eigenvalues of GβE at the edge for β ∈ {1, 2, 4} [104,105] we obtain that the mean density for the fermion model takes the scaling form in the large N limit where the scaling functions are with κ = 2 2/3 . The (smooth) linear statistics for general β was studied in [106] in the limit towards the bulk. For the FCS for general β however explicit formulae are still lacking in the edge region.
Concerning the variance of the particle number, the scaling form (58) can be extended to any β, see [41,42] However the explicit form for general β is unknown at present. The asymptotic behavior in the limit towards the bulk,â → −∞ can also be extracted from the results in [77], see Appendix F for details. This leads to the asymptotic behaviors for β ∈ {1, 2, 4} which matches with the bulk result (40) that we obtained above. The leading order (logarithmic) term in (65) was conjectured in [41] for any β based on the expected matching with the bulk. From our conjecture in Section 2.4 we can now predict that (65) holds for general β, with the constant c β given in (17).
Concerning the cumulants of the particle number of order three and higher, these have been obtained in the limit towards the bulk (â → −∞) for β ∈ {1, 2, 4} in Ref. [77]. In Appendix E we have verified that our prediction for the higher cumulants (20) for general β match perfectly with the results from Ref. [77] for β ∈ {1, 2, 4}. This provides a non trivial check of our conjecture, and involves not so well known identities among Barnes functions.
The above results were obtained for fermions in the harmonic potential, associated to the GβE. In RMT it is known that there is a universality at the soft edge, hence we expect the same behavior to hold for the fermion model associated to the WLβE (line 3 in Table 1) at its soft edge (i.e., near x + e = √ 2βN ). On the other hand, for that model near x − e = 0 (and for the Jacobi box) the universality is called the hard-edge and described by the Bessel stochastic operator [107].
As in the case β = 2 it is reasonable to conjecture that for general β the above results extend to any smooth confining potential, e.g. V (x) ∼ |x| p with p > 0, near the edge, the functions V β being thus universal. The simple picture is that V (x) near the edge can be approximated by a linear potential in all cases. One can also expect that the higher cumulants will also take a universal scaling form, being a non-trivial function ofâ.
Note that the effect of short range interactions at the edge was discussed in Ref. [108] and found to be subdominant within the model studied there. It was also noticed there that the case of 1/x 2 interactions lead to a new universality class for general β.

Bosonisation and Luttinger liquid
It is well known that interacting fermions in one dimension can be described by the effective theory of the Luttinger liquid (LL) based on the bosonisation method, for a review see e.g. [109]. It provides a hydrodynamic description which in its simplest form is valid in the absence of an external potential. At equilibrium, and for spinless fermions, it depends only on two parameters, the mean density ρ 0 and the dimensionless Luttinger parameter K, with K = 1 for noninteracting fermions, K < 1 for repulsive interactions and K > 1 for attractive interactions. The dynamics also depends on the "sound velocity" v F (Fermi velocity for free fermions). This is based on the description in terms of the phase field ϕ(x) defined such that ρ(x) = − 1 π ∂ x ϕ(x), which at large scale is described by a Gaussian theory. This theory can be extended in the presence of a potential which varies very slowly on scales of the order of the inter-particle distance 1/ρ(x). Many recent studies have addressed the case of inhomogeneous bosonisation where v F , K and ρ 0 may become slowly varying functions of the position x [51][52][53] (see also [110]).
Let us recall that for a LL the density correlations are given by while the correlation function of the fermionic field (which is the analogue of the kernel in the case of noninteracting fermions) reads These formulae (66) and (67) are valid for ρ 0 x 1. Here A 1 in (66) and C 0 in (67) represent the leading behaviors at large ρ 0 x, while the terms A m , m ≥ 2 and C m , m ≥ 1 represent the contributions of higher harmonics (often neglected in LL studies). For noninteracting (free) fermions, K = 1, C 0 = 1 π and all C m = 0 for m ≥ 1, and the expression in (67) becomes exact. In this case, this is precisely the sine kernel sin(πρ 0 x)/(πx). In the presence of interactions we see that the correlation function of the fermionic field in (67) in the ground state now decays at large x as with a non-universal prefactor. Consider now the model of interacting fermions on the circle (second line in Table 1). One can predict that it corresponds to a Luttinger liquid with parameter K = 2/β. Indeed the density correlations were calculated for the CβE in Ref. [111] and one can check that formula 4.11 there agrees with the prediction of the LL theory (66), with the choice K = 2/β, up to subleading terms (which for each harmonic decreases faster by a factor 1/|x|). As mentioned in proposition 13.2.4, p. 604 of [68] (see also [112,113]) this asymptotics was established for β an even integer. However the value of K can be inferred already from the second term in (66), i.e., from the coefficient of the long range decay ∼ 1/x 2 , which in fact can be obtained by electrostatic arguments and linear response theory from the Coulomb gas representation [84,114]. Note that the identification K = 2/β was also noted in [108] (see Appendix there) and in Ref. [115] where an approximate formula for more general powerlaw interactions was also obtained. We thus expect that all the universal properties of the Luttinger liquid with this value of the parameter will hold for the model in second line in Table 1 for fermions on the circle. For instance the variance of the number of fermions, N [a,b] = 1 π (ϕ(a) − ϕ(b)), can be computed from the correlator of the phase field given in [109] (see also [116]) as Fermions' domain where |a − b| 1/k F , where in the general interacting case k F is defined as πρ 0 and is unrelated to v F (while k F = v F /m in the noninteracting case). Substituting K = 2/β in (69) this agrees with the leading logarithmic term in our prediction (19), although it is not accurate enough to predict the O(1) term.
In the presence of both interactions and an external potential, an inhomogeneous bosonization approach was recently developed which aims to calculate correlations in the bulk using conformal field theory methods [51][52][53]. The potential induces a spatial dependence of the LL parameters, K(x), v F (x) and ρ 0 (x). In the case of free fermions this approach can be compared with the exact results of Ref. [49] that we presented in Section 2. 1. The correspondence appears to be, from Eq. (20) of [51] y) is defined there and θ(x) was defined in (25). It would be interesting to extend these results to the present models with K = 2/β. Since it is natural to assume, because of the long range nature of the interactions, that K = 2/β is independent of x, a similar description will hold for more general potentials than the ones considered here.

More general models
So far we have focused on the models in Table 1, with ground state wave functions of the form (5) involving one and two-body factors and are related to random matrix models of the form (6). There is in fact a larger class of models for which the ground state wave-functions are still of the form (5). The study of such models was initiated by Sutherland and Calogero [65,67,117] and extended in Refs. [118][119][120]. We recall how models with this property are constructed using a slightly more general approach, and give a list of corresponding Hamiltonians in the Appendix A. In addition we show that some of these models are related to other interesting random matrix models. Some of them were studied in the RMT literature. The models presented in this section are summarized in Table 2.
Fermions on the circle. The first interesting extension corresponds to fermions on the circle, i.e x ∈ [0, 2π] with periodic boundary conditions, in the presence of an external periodic potential. It generalizes the models of the second line of the Table 1 which are related to CβE (in the absence of potential). It corresponds to the Hamiltonian (4) with the external potential V (x) and two-body interaction W (x, y) given by W (x, y) = 1 16 The ground state wave function is of the form (5) for any N , with v(x) = g cos x and w(x, y) = −β log | sin x−y 2 | and the ground state energy E 0 is given in (116). The quantum probability is (up to a normalization) Interestingly, for β = 2 this probability is identical to the one studied at large N by Gross and Witten and independently by Wadia in the context of lattice gauge theories [78,79], and later on in combinatorics [121]. The weight (72) corresponds to the joint PDF of the eigenvalues λ j = e ix j of a matrix model with a probability measure on the N × N unitary matrices ∝ e g 2 Tr(U +U † ) dU . Remarkably, in the present context this corresponds to noninteracting fermions. The mapping is summarized in the first line of Table 2.
As we now show, this mapping allows us to shortcut the Coulomb gas method used in Refs. [78,79] and obtain the eigenvalue density in a simpler way. For β = 2 the external potential for the fermions reads where we have definedg = g/N which is the important parameter that is kept of O(1) in the large N limit. In Fig. 5 we show a plot of V (x) for various values ofg. Since there is no interaction W = 0 for β = 2 in the large N limit the fermion density can be obtained from the LDA formulaρ where µ is the Fermi energy and is determined by the normalization condition 2π 0 dxρ(x) = 1. From (73) and Fig. 5 we see that there are two cases depending ong. Forg < 1 there is a single maximum V max of the potential V (x) which is attained for u = cos x = 1, i.e., x = 0, with V max = N 2 8 (2g −g 2 ). Forg > 1 there are two degenerate maxima of the potential V (x) which are attained for u = cos x = 1/g, of values V max = N 2 8 . As we now discuss this change of behavior of V (x) results in two distinct phases: forg < 1 (weak coupling) the Fermi energy is above the maximum of the potential and the density is everywhere positive. Forg > 1 (strong coupling) the Fermi energy is below the maximum and the density has a restricted support.
Since the Fermi energy µ is itself determined by the normalization condition this is a non trivial transition. In the weak coupling phase one finds, as we show below, that For that particular value, one sees that the expression inside the square root in (74) becomes a perfect square leading toρ which is automatically normalized to unity, with a support x ∈ [0, 2π]. Hence (75) is the correct value for µ. Since the density must be positive, this solution is acceptable only for g < 1.
Forg > 1 we note that the potential V (x) has a local minimum at x = 0 of value V min = N 2 8 (2g −g 2 ). In this strong coupling phase, one finds As can be seen on the Fig. 5 this value of µ is such that the system has a single support in all phases. For this value of µ we obtain from (74) which has now a restricted support [x − , x + ] where the edges are x − = 2arccos 1 g and x + = 2π − x − . One can check the normalization which shows that (77) is the correct value of µ. Forg → 1 + one has x − 2 √g − 1, and for g → +∞ one has x + → π − (all fermions are around x = π). Forg = 1 the formulae (76) and (78) become identical.
The phase transition in the density in the above formula recovers the results obtained in [78,79] by a different method, upon the identificationg = 2/λ (and x → x + π) from the notations of [78]. In these papers the partition function (i.e., the normalization amplitude of the probability measure in (72)) was computed and shown to exhibit a third order phase transition atg = 1 (for a recent review, see [103]). In the fermion system this transition at g = 1 can be seen as a freezing transition for the Fermi energy as a function of the coupling strengthg (see Figs. 5 and 6) µ = N 2 8 ,g < 1 (weak coupling) N 2 8 2g −g 2 ,g > 1 (strong coupling) (80) This transition coincides with the opening of a gap in the bulk in the fermion density. A similar transition has been recently studied by us for fermions in an inverted parabolic potential [122], and the correlation kernel at the transition was explicitly computed. However, in the present situation the critical behavior is expected to be different. Indeed the shape of the potential at criticality around Finally note that the model for β = 2, i.e., for interacting fermions in (70) is also of interest. Using CG arguments, its densityρ β,g (x) is obtained from the density for β = 2 by a simple rescaling, i.e.,ρ β,g (x) =ρ 2, 2 βg (x) found above. The transition then occurs forg = β/2. The correlation functions are expected, however, to depend on β and remain to be explored.
Fermions on the half-line. The second interesting extension generalizes the models on the third line of Table 1 which are related to the WLβE. The interaction potential W is the same as in Table 1, but the potential V is more general The ground state wavefunction is of the form (5) with v(x) = c 0 x 2 + c 1 x 4 + c 2 log x. It corresponds to a matrix model of the form (7) upon the map λ(x) = 2 β x 2 with a matrix potential The mapping is summarized in the second line of Table 2.
Fermions in a box. The next extension generalizes the models on the fourth line of the Table 1 which are related to the JβE. The interaction potential W (x, y) is the same as in Table  1, but the potential V (x) contains additional cos x and cos 2x terms, see formula (119). The ground state wave function has the form (5) with v(x) = c 1 log sin x 2 + c 2 log cos x 2 + c 3 cos x. It corresponds to a matrix model of the form (7) upon the map λ(x) = 1 2 (1 − cos x) and matrix potential The mapping for this model is summarized in the third line of Table 2. For c 3 = 0 it recovers the JβE. For c 1 = −1 and c 3 = 0, this matrix model was studied in [150][151][152] and its density was calculated using the Coulomb gas method. We show in Appendix B that this result agrees with the LDA. Finally there are some models not related to Table 1. Hyperbolic models. The simplest example are fermions on the real line with the two-body interaction potential Its ground state wave function is of the form (5) with a two-body term w(x, y) = −β log | sinh 1 2 (x− y)|. For normalizability of (5) one needs a confining potential. The most general family consistent with (84) is The mapping for this model is summarized in the fourth line of Table 2. In the case c 2 = 0 of the Morse potential this relation to the Wishart model was also obtained in [123]. For c 2 = 0 the calculation of the mean density σ(λ) was performed using Coulomb gas methods for some values of the parameters in [124]. We show that this result agrees with the LDA in Appendix B. Note that this matrix model was also studied in various contexts in Refs. [125,126].
Another interesting example in this class of hyperbolic models corresponds to a ground state wave function of the form form (5) with v(x) = ax 2 and w(x, y) = −β log | sinh 1 2 (x − y)|. In that case there is an additional repulsive two-body interaction δW on top of the interaction (84), of the form δW (x, y) ∝ −aβ (x − y) coth 1 2 (x − y), which never vanishes for any value of β (the model is always interacting). This model is related to the Stieltjes-Wigert β-ensemble (SWβE) [119,[127][128][129][130] which was studied in the context of Chern-Simons theory in high energy physics [131] and of non-intersecting Brownian bridges [130,132,133]. The correspondence is through the map (see the discussion below Eq. (125) and the matrix potential reads whereã = 2a/β. The mapping for this model is summarized in the sixth line of Table 2. For this matrix model the joint PDF of the eigenvalues is determinantal for β = 2, since the model becomes bi-orthogonal [134][135][136]. In the limit of large N , scaling a = O(N ), the eigenvalue density σ(λ) is known [130,131,133] where u = N/(2ã) = O(1). From the CG arguments this density is in fact independent of β.
Hence we obtain the fermion density for the associated quantum model for any β as ρ(x) = N e u/2 e x σ(e u/2 e x ).
Finally, there are two more hyperbolic models for fermions, one on the positive half axis which corresponds to the fifth line in Table 2, and the second one, which maps to the Cauchy random matrix ensemble, and corresponds to the seventh line in Table 2. These models are described in the Appendix A.
Let us close this section by indicating yet another family of quantum models where the interaction W (x, y) is a sum of a harmonic attraction ∝ (x − y) 2 and of the inverse square interaction β(β−2) 4(x−y) 2 . The first case is in an external potential V (x) ∼ x 2 . The second is related to a quartic matrix model V 0 (λ) = c 2 λ 2 + c 4 λ 4 and corresponds to a fermion model with a polynomial potential with terms x 2 , x 4 , x 6 . These models are described in the Appendix A.
To relate to the main focus of the paper, i.e., the counting statistics, let us point out that many models presented in this Section are noninteracting for β = 2. In that case, the methods of [49] summarized in the Section 2.1 can be applied to obtain the variance of the number of fermions in an interval. Upon scaling properly the parameters of the model with β one can relate the variance of the interacting model to the one for β = 2 by similar relations as in (18), with the same constants (17). Our conjecture for the higher cumulants should also apply.

Discussion and conclusion
In summary, we calculated the counting statistics for several models of N 1 interacting spinless fermions in their ground state in one dimension confined by an external potential, see Tables 1 and 2. The interactions are of the general Calogero-Sutherland type, and depend on the parameter β, where β = 2 corresponds to the noninteracting case. We have emphasized the connections to random matrix ensembles, where β is the Dyson index. We found that the variance of the number of fermions in a macroscopic interval [a, b] in the bulk of the Fermi gas grows with N as A β log N + B β + o(1). We obtained explicit formulae for A β and B β , which depend on a, b, on the type of interaction and on the shape of the confining potential. These results were obtained by explicit calculations for β ∈ {1, 2, 4} and from a conjecture that we formulated for general β. This conjecture extends to the higher-order cumulants of the distribution of N [a,b] . They are O(1) at N 1 and are predicted here to be given by (51). Remarkably, this result is universal: it does not depend on the confining potential. This is because the conjecture states that the short scales determine the O(1) part of the fluctuations of the particle number. We have obtained a few "smoking gun" tests for this conjecture. First we have shown that it matches in a highly nontrivial way, near the edge of the Fermi gas, with recent results from the mathematics literature [77]. Second we have shown that our analytical predictions are in very good agreement with our numerical simulations. In addition we have shown that the leading term A β is in agreement with the predictions from the Luttinger liquid theory with parameter K = 2/β.
Finally, we formulated a general approach for obtaining mappings between interacting fermion models in one dimension in their ground state and random matrix models (or, more generally, models of classical interacting particles confined by an external potential in thermal equilibrium). We applied this approach and found several such mappings. In particular we found a surprising mapping of the famous RMT Gross-Witten-Wadia model from high energy physics onto noninteracting fermions in an external potential on a circle. The simple application of the LDA allows to obtain the mean fermion density in that case, and recovers results known for this model obtained by more involved Coulomb gas methods. In turn, we have shown that these Coulomb gas methods can be used to study interacting fermions in a trapping potential. We exploited these mappings to obtain the mean fermion density for these models for general interaction parameter β by relating them to the noninteracting case β = 2. Similarly, we argue that the counting statistics in these models can be calculated by relating them to the noninteracting case.
Our results hold also for Dyson indices 0 < β < 1 which, although meaningless in the fermion systems on which we focused here (since for fermions β ≥ 1), are meaningful for the RMT ensembles. The scaling limit β ∼ 1/N has generated much interest recently, and it would be interesting to study the counting statistics in this limit [137][138][139][140][141].
Among the connections unveiled in this paper, e.g., with the models in Table 2, many interesting questions remain to be explored. In particular one may wonder whether the universality of the higher cumulants of the fermion number can be extended to more general interacting models, and whether one can derive formula for the variance in more general potentials. In particular, it would be interesting to test this universality when perturbing the interaction term away from the Calogero-Sutherland type studied here.
For noninteracting fermions, the counting statistics is connected to the bipartite entanglement entropy (EE) of the subsystem D with its complement D [12][13][14][15]. Given the results of the present work, it would be interesting to search for similar (perhaps approximate) connections for interacting fermions in order to calculate the EE.
Finally, it would be interesting to extend our approach to higher dimensions. In particular, there is a known mapping between noninteracting fermions in a 2d rotating harmonic trap and random matrices of the complex Ginibre ensemble [16,17]. It remains a challenge to extend this mapping to more general cases.
Acknowledgements PLD thanks Y. V. Fyodorov for an earlier collaboration on related topics. We thank A. Borodin and P. J. Forrester for useful correspondence. We thank D. S. Dean and C. Salomon for interesting discussions. We thank T. Bothner for useful comments on the manuscript. We thank M. Beau for pointing out the recent references [142,143] about ground-states in Calogerotype models and their extensions in higher dimensions. NRS acknowledges support from the Yad Hanadiv fund (Rothschild fellowship). This research was supported by ANR grant ANR-17-CE30-0027-01 RaMaTraF.
A Interacting fermion models with ground state of the form (5) and mappings to RMT In this Appendix we recall the construction of quantum Hamiltonians with two-body interactions in one dimension (4), whose ground state wave function has itself a two-body form as in Eq. (5). This question was pionneered by Calogero [67] (following Sutherland [65,117]) and extended in Refs. [118][119][120]142,143]. In some cases these models are also fully integrable (i.e., their full eigenspectrum is known), see e.g. [66,68,144,145]. Here we also discuss the construction of the ground state in the light of the connections to random matrix ensembles. In particular we perform a search for models using the map λ(x) which relates RMT to fermions.

A.1 Schrödinger equation and general conditions for two-body-only interaction
Consider the following unnormalized wave function, Ψ 0 ( x) = e −U ( x)/2 (defined up to a sign in an ordered sector), where U ( x) = i v(x i ) + i<j w(x i , x j ) has the two-body form (5). A necessary condition for it to be the ground state of the two-body Hamiltonian H N in (4) with energy E 0 is that H N Ψ 0 ( x) = E 0 Ψ 0 ( x). Substituting and multiplying by e U ( x)/2 on both sides one gets where T n denotes a term which is naively n body. Here, and in the following, we use the We recall that w(x, y) is a symmetric function and denote by subscripts the order of its partial derivatives. These terms are where we have splitted the term 1 8 i j =i w 10 (x i , x j ) k =i w 10 (x i , x k ) into the term j = k (in T 2 ) and j = k (in T 3 ). In the cases that we will study, these terms will drastically simplify and turn out to be constants (and sometimes zero). The ground state energy E 0 will be determined, as a result.
To obtain a two-body Hamiltonian we must thus impose the condition that the three-body interactions are absent. This amounts to a condition on w(x, y) so that T 3 can be written as two-body term, or a one-body or a constant. We will search for solutions to this condition in two possible forms w(x, y) = w(x − y) and w(x, y) = −β log |λ(x) − λ(y)|. Asking that T 3 is a constant, or one-body, then allows for a systematic search. This leads to a set of quantum model with two-body interactions W (x, y) = (W (1) + W (2) )| 2body , with a specific family of interactions T 2 | 2body which vanish for β = 2, while T 2 | 2body depends on v(x). For some specific choices of v(x) which we identify T 2 | 2body = 0, which lead to simpler quantum models in an external potential which become noninteracting for β = 2.
In the next section we make the list of the models which are obtained by this method, and in the following section we explain how one searches for these models.
Remark: The term T 1 has the form of potentials from supersymmetric quantum mechanics.

More generally the above equation (90) is equivalent to
For repulsive interactions β > 2, the mappings described here are expected to hold for bosons too. For bosons, it is the repulsive interaction that causes the many-body wave function Ψ 0 to vanish at x i = x j for i = j.

A.2 Families of models
We consider here different kinds of models. Some are defined on the real line (or the half-line) and require a confining potential v(x) in order for Ψ 0 to be normalized. The others are called "periodic" models, and defined either on the circle or an interval, in which case v(x) may be chosen to be zero. We recall that when there is a mapping x → λ(x) between the fermions models with potential v(x) and a matrix model (6) with a matrix potential V 0 (λ), the relation between the two potentials reads Note that w and v in (5) are defined up to an irrelevant additive constant which can be absorbed into the normalisation of Ψ 0 . Depending on v(x), one may also extract a one-body part from W (x, y) and add it to V (x), and extract constant parts from W, V and add them to −E 0 , where E 0 below denotes the ground state energy.
Logarithmic models. In this class, the first set of models is, for x, y on the real axis w(x, y) = −β log |x − y| , λ(x) = x , T 3 = 0 (100) In this set of models the only normalizable choice of v(x) which does not contribute to the two-body interaction W is v(x) = ax 2 . It corresponds to the quantum model [146] V (x) = a 2 2 This corresponds to the GβE, for which the canonical choice, given in the text, is a = 1, λ(x) = 2 β x and V 0 (λ) = βλ 2 /2. For the noninteracting case β = 2 (setting a = 1) one recovers that E 0 is the sum of the energies of the single-particle states up to the Fermi energy, The second set of models is, for x, y on the positive real axis [154] w(x, y) = −β log |x 2 − y 2 | , λ(x) = x 2 , T 3 = 0 (103) In this set of models the only normalizable choice of v(x) which does not contribute to the two-body interaction is v(x) = c 0 x 2 + c 1 x 4 + c 2 log x which corresponds to the quantum model [146] V (x) = 2c 2 1 x 6 + 2c 0 c 1 x 4 + This contains the case of the WLβE with the canonical choice, given in the text which leads to For β = 2 one recovers the energy E 0 = N −1 n=0 (2n + 1 + γ). However there is a larger class of potentials which correspond to matrix models with matrix potentials V 0 (λ) = c 0 Periodic models. In this class, the first set of models is defined on the circle with px ∈ [0, 2π[ It contains the CβE which is obtained for v(x) = 0. The canonical choice given in the text is p = 1. One can check that for β = 2, the ground state energy is exactly equal to the sum of the energies of the single-particle states, e.g. for N odd one has E 0 = 2 24 .
In this set the only choice of v(x) which does not generate a two-body interaction is v(x) = b cos(px) (up to translations on the circle), which leads to the quantum model on the circle W (x, y) = p 2 16 For β = 2 this is the Gross-Witten-Wadia model discussed in the text.
Hyperbolic models. In this class, the first set of models is defined on the real axis It is thus equivalent to a matrix model with the matrix potential V 0 (λ) such thatṽ(x) = V 0 (e px ) − px.
For the choice v(x) = ax 2 this model is related to the Stieltjes-Wigert β ensemble [119,130] as discussed in the text, where we have also used the parametrization (w,ṽ) to obtain Eq. (87).
In this set the only choice of v(x) which does not contribute to the two-body interaction is v(x) = c 0 x + c 1 e px + c 2 e −px . This leads to the quantum model [146] V (x) = 1 8 which has a potential of the (generalized) Morse type. Note that the ground state energy behaves as ∝ −N 3 , which is surprising for a repulsive interaction (say β > 2). This is because the confining potential (necessary for normalization) has a deep minimum at energy ∝ −N 2 . Note that it corresponds to a matrix model with matrix potential The second set of models is defined on the positive half line and corresponds to the choice w(x, y) = −β log | cosh px − cosh py|, which is equivalent to the choice In this set the only choice of v(x) which does not contribute to a two-body interaction is v(x) = c 1 log sinh px 2 + c 2 log cosh px 2 + c 3 cosh px, which leads to the quantum model [146] V It corresponds to a matrix model with matrix potential The mapping for this model is summarized in the fifth line of Table 2.
Finally there is a third set of models, defined on the line or on the positive half-line and which corresponds to the choice w(x, y) = −β log | sinh px − sinh py|, which is equivalent to the choice In this set the only choice of v(x) which does not contribute to the two-body interaction is v(x) = c 1 arctan sinh px + c 2 log cosh px + c 3 sinh px which leads to the quantum model [146] V (x) = c 2 1 − c 2 (c 2 + 2) Note however that for normalizability on the whole axis one needs c 3 = 0. For c 3 = 0 the potential V (x) is known as the hyperbolic Scarf potential [155]. It corresponds to a matrix model The mapping for this model is summarized in the seventh line of Table 2. In the case c 3 = 0 this model is the generalized Cauchy beta ensemble (CaβE) studied, e.g., in [156,157]. The joint PDF of eigenvalues has the form with a = 1+c 2 2 and b = −c 1 /2. In the case c 1 = 0 (i.e., b = 0) and a = β 2 (N − 1) + 1, the CaβE is related to the circular ensemble CβE via the stereographic projection e iθ = 1+iλ 1−iλ . As a result the exact density σ N (λ) is known for any finite N from the fact that it is uniform on the circle [68,157] and it is given by independently of N and β. This mapping does not relate however the Schrödinger operators on the circle and the line, so it does not extend to the quantum model. In fact the potential in the fermion model associated to the case c 1 = c 3 = 0 is the Pöschl-Teller potential, as can be seen from (142) (setting p = 1) For β = 2 (noninteracting fermions) the LDA approximation at large N for the fermion density, is ρ(x) = 1 π 2µ + c 2 (c 2 +2) 4 cosh 2 x . It is easy to check that this formula is compatible with the exact result , which holds for c 2 2N at large N . This determines the value of the Fermi energy as µ N 2 0, which means that in the ground state, the potential well which has only a finite number of energy levels, is almost full.
Elliptic models. In all the above models T 3 was a constant. There is a more general family of models for which T 3 is a sum of one-body terms. They are such that w(x, y) = −β log |λ(x) − λ(y)| and λ(x) solution of In that case λ (x) , which leads to a potential V (x) = 1 8 (N − 1)(N − 2)u(x). For C = 0 one recovers the models discussed above (and u(x) is then a constant). For the general case C = 0 the solutions of (149) are of the form λ(x) = a + b P(x; g 2 , g 3 ) where a = −A/C, b = 12/C and g 2 = − C 6 (B − A 2 2C ) and g 3 is an arbitrary constant. Here P is the Weierstrass function [158] (see also next Section). This leads to V (x) = β 2 2 (N − 1)(N − 2)P(x; g 2 , g 3 ). The two-body term T 2 then leads to the interaction W (1) which appears to be quite complicated. We will not further study these solutions here.
Models with quadratic interactions. A last example is the following fermion model defined on the real axis For general a, one has that T 3 is the sum of two-body terms. There are two interesting special cases. The first one is the case where v(x) = c 2 x 2 which leads to [146] W (x, y) = a 2 + a 2 2 Another special case is a = 0, for which T 3 is actually a constant (this model also belongs to the first class in Eq. (100)). In this case, one has w(x, y) = −β log |x − y| and v(x) = c 2 x 2 + c 4 x 4 which leads to [146] W (x, y) = V (x) = 2c 2 4 x 6 + 2c 2 c 4 x 4 + This quantum model is thus related via λ(x) = x ∈ R, to the random matrix model with a quartic matrix potential V 0 (λ) = c 2 λ 2 + c 4 λ 4 .
This model is well known in RMT [147,148,159,160]. The mapping for this model is summarized in the eighth line of Table 2.

A.3 Search for models
Let us summarize how one searches for models. One imposes the condition that the threebody interactions are absent, i.e., that T 3 can be written as a two-body term, or a one-body term or a constant. We first consider the case when T 3 in Eq. (97) is simply a constant. It means that for all x, y, z t 3 (x, y, z) := w 10 (x, y)w 10 (x, z) + w 10 (y, x)w 10 (y, z) + w 10 (z, x)w 10 (z, y) = ±q 2 , in which case Note that this does not involve v(x) hence for now v(x) is arbitrary.
A first series of model is obtained by considering w(x, y) = w(x−y) with w even. Inserting into (160) and taking z, y → x one sees that there is no differentiable solution at x = 0, i.e., with w (0) = 0 since w (x) is an odd function. Hence one needs w (x) to diverge at x = 0. One thus writes w (z) = 1/g(z) with g(z) an odd function. Setting z = x + , one must have, regrouping the terms We see that the only possibility (for a smooth g(x) at generic non zero x) is g( ) = O( ), i.e., a simple pole for w (x), and taking → 0 one finds The solutions are g(x) = 1 q tan(qxg (0)), g(x) = 1 q tanh(qxg (0)), and g(x) = g (0)x. Defining g (0) = −1/β and q = − β 2 p one obtains w(x, y) = −β log | sin( p 2 (x − y))| for the (+) branch, w(x, y) = −β log | sinh( p 2 (x − y))| for the (−) branch, and w(x, y) = −β log |x − y| for q = 0. Here β and p are for now an arbitrary parameters. One checks that indeed they satisfy the condition (160) (providing not so trivial trigonometric identities). These are the models respectively (110), (125) and (100). For these models the two-body term from T 2 leads to the interaction potential W (1) In the presence of a potential v(x) there is generically another two-body term from T 2 in Eq. (96). It reads W (2) (x, y) = 1 4 (v (x) − v (y))w (x − y) which leads to the interaction potentials in (110), (125) and (100). For each of these three models there is a unique family of exceptional potentials v(x) for which W (2) (x, y) is either one-body or constant. To search for them one imposes the necessary condition ∂ x ∂ y W (2) (x, y) = 0 and solves it for y → x. This leads to the models (102), (114) and (128).
Interestingly, the models found so far are closely related, up to some change of variable λ = λ(x), to the logarithmic interaction which appears naturally in the β random matrix models of the form (6), with an a priori arbitrary matrix potential V 0 (λ). Thus it is natural to search for w(x, y) parameterized in the form w(x, y) = −β log |λ(x)−λ(y)|, where again, for now, v(x) is arbitrary. Note that there is some redundancy, since e.g. the problem with λ(x) = 1/λ(x) is equivalent to w(x, y) = −β log |λ(x) −λ(y)| and v(x) → v(x) + (N − 1)β log |λ(x)|. One now imposes to satisfy (160) with w 10 (x, y) = −β λ (x) λ(x)−λ(y) and we restrict here to λ(x) real. Substituting and taking successively the limits z → x and y → x, one arrives at the necessary condition β 2 λ (x) λ (x) = ±q 2 . It is convenient to parameterize these models using the parameter p defined now via q 2 = β 2 p 2 [149]. The general solutions are λ(x) = ax + bx 2 (for q = 0), λ(x) = a 1 cos px + a 2 sin px (for the − branch) and λ(x) = a 1 cosh px + a 2 sinh px for the (for the + branch). The case λ(x) = x recovers (100), while λ(x) = x 2 gives the new family (103) (and upon a translation in x one can always reduce to one of these cases).
Similarly for the periodic model one can always choose a 2 = 0 by translation, which leads to (117). Finally for the hyperbolic solutions one can always reduce to either a 1 = a 2 , which leads again to (125), or to a 2 = 0 (whenever a 1 > a 2 ), which leads to (134), or to a 1 = 0 (whenever a 1 < a 2 ), which leads to (140). One can check that T 3 is indeed a constant for each of these models, again via non trivial trigonometric identities. With the chosen parameterization w(x, y) = −β log |λ(x) − λ(y)| the term T 2 leads to an interaction λ(x)−λ(y) ), which simplifies for these models as given in (100), (103), (117), (125), (134), (140).
Next we turn to the condition that the three-body term in (160) is the sum of one-body terms, i.e., t 3 (x, y, z) = u(x) + u(y) + u(z). Substituting and taking successively the limits z → x and y → x one arrives at the necessary condition β 2 λ (x) λ (x) = 3u(x). Next we expand in y − x and to second order we obtain another necessary condition = Cλ (x). Multiplying by λ (x) and integrating once more leads to the simple condition λ (x) = B + Aλ(x) + C 2 λ(x) 2 , i.e., the relation given above in Eq. (149), where A, B and C are integration constants.
We have not explored in full generality the condition that T 3 is a two-body term, i.e., t 3 (x, y, z) defined in (160) be a sum of two and one-body potentials. This condition leads to the nonlinear, nonlocal partial differential equation ∂ x ∂ y ∂ z t 3 (x, y, z) = w 21 (x, y)w 11 (x, z) + w 11 (x, y)w 21 (x, z) + w 21 (y, x)w 11 (y, z) +w 11 (y, x)w 21 (y, z) + w 21 (z, x)w 11 (z, y) + w 11 (z, x)w 21 (z, y) = 0 , (167) whose study we leave for future research. In the case of the form w(x, y) = w(x − y) this was done by Calogero [67] who found that the general solution must obey w (z) = a P(z, g 2 , g 3 ) where P is the Weierstrass P function, i.e., a solution of P = 6P 2 − g 2 2 and (P ) 2 = 4P 3 − g 2 P − g 3 [158]. Integrating twice the general solution is thus in that case, w(z) = −β log(σ(z; g 2 , g 3 )) + bz 2 where σ(z; g 2 , g 3 ) is the σ-Weierstrass function. Thanks to non trivial identities involving Weierstrass functions [67], the absence of three-body term is indeed obeyed. Note that it is natural to set a = β since σ z at small z, hence w(z) is again of the logarithmic type at small z. The resulting quantum interaction W (x, y) can be expressed as a polynomial in terms of Weierstrass functions [67].
Here we only give the simple example (150). In fact, the only solution with w(x, y) = w(x − y) and w (0) finite is the quadratic form w(x, y) = a(x − y) 2 .

B Mean fermion density: LDA versus Coulomb gas
In the absence of interactions, β = 2, the LDA prediction for the mean fermion density in the large N limit is (with unit normalization) On the other hand, when there is a map λ(x) which maps the ground state of this model to the joint PDF of the eigenvalues of a RMT ensemble of the form i.e., when the matrix potential in (6) is scaled as V 0 (λ) = β 2 NṼ 0 (λ), then it is possible to use the CG method to obtain the eigenvalue density σ(λ). It is given as the optimal density which minimizes the CG energy functional under the constraint that dλσ(λ) = 1, which we (abusively) also denote σ(λ). The connection between the two isρ(x) = λ (x)σ(λ(x)). Since the CG density is independent of β this allows to obtain the fermion density for any β.
We have discussed this connection in the text on the example of the Gross-Witten-Wadia model. Here we give some more details for the other cases. The computationally difficult part is to determine the Fermi energy µ.
Hyperbolic model. Consider the model (84) discussed in the text with potential V (x) given in (85), which has three parameters c 0 , c 1 , c 2 . We will scale them as c j = N β 2c j . From (86), it corresponds, in the large N limit, to the matrix potential (dropping subleading terms at large N )Ṽ 0 (λ) =c 1 λ +c 2 λ −1 + (1 +c 0 ) log λ. In [124] the minimization equation was solved in the casec 1 = 1, 1 +c 0 = −1 and it was found that In [124] the parameters a, b, c are given as a function of µ 1 =c 2 . On the other hand the LDA prediction for β = 2 with λ(x) = e x gives in the more general case of the potential (85) Let us recover (171) from this result. Pluggingc 1 = 1,c 0 = −2 into (172) yields Now, requiring that the expression under the square root in (173) can be written in the form for some a, b, c, one obtains by comparing the coefficients of powers of λ on both sides of the equation, the following relations 6 = a + b − 2c, 8μ = −ab + 2ac + 2bc − c 2 , −2c 2 = −2abc + ac 2 + bc 2 ,c 2 2 = abc 2 , (175) whose solution, in terms of v = √ ab and u = a/b, is v = 2u 3u 2 − 2u + 3 in agreement with [124], and the (rescaled) fermi energỹ which can be expressed in terms ofc 2 alone using the equations in (176). In particular, the agreement with [124] ensures that σ(λ) is correctly normalized, which shows that the form (174) is indeed correct.
Fermions in a box. This is the case of the matrix potential (83). Here λ(x) = 1 2 (1 − cos x). The LDA prediction is given by (169) where V (x) is given in (119) (with β = 2) and µ is determined through the normalization ρ(x)dx = 1. In general, this expression forρ and the calculation of µ are very cumbersome. However, choosing p = 1 and β = 2 in (119) and assuming c 1 = O (1) , c 2 = O (1) , c 3 = O (N ), the potential is simply to leading order for large N . We note that up to the additive constant c 2 3 16 this potential coincides (at N 1) with the potential (73) with 2N particles, if one identifiesg ↔ c 3 2N . Therefore, the density predicted by the LDA can be immediately deduced from Eqs. (75)- (78) with 2N particles, by adding a factors of 2 toρ due to the difference between the domains [−π, π] for Gross-Witten vs. [0, π] in the present case. The result is: and the associated fermi energies are It is assumed above that c 3 > 0, but flipping the sign of c 3 is equivalent to transforming This can be compared with the predictions of [150][151][152] who studied the model (6) with matrix potential (83) with (in our notations) They obtained the density which, for β = 2, leads to in agreement with (179).

C Number variance for the harmonic trap
Here we calculate the number variance for interacting fermions described by the model (3), which corresponds to random matrices in the GβE, thereby obtaining Eqs. (40) and (14) of the main text. Let us first consider Var N [a,∞) . We aim to calculate the double integral (31) in the large-N limit, by approximating C(x, y) 0 if either x or y are not in the bulk, and for x and y both in the bulk, plugging in C (x, y) from (32) for x near y, and (38) for x far from y (this procedure works because there is a joint regime of validity for both of these approximate expressions for C(x, y)). Thus Var N [a,∞) I 1 + I 2 where where we have chosen some cutoff ξ such that 1 √ N ξ 1, which justifies the approximation ρ N (x) ρ N (a) that we made in the integral I 2 .

D Number variance for the WLβE and the JβE
In this Appendix we present a detailed derivation of Eq. (48) and we give the result for the variance of the fermion number for the models in the third and fourth line of the Table 1, which are not already given in the text.
In [49] we found the number variance for the WLβE with β = 2 for a semi-infinite interval: where µ = 2N + γ + 1 and λ 2 = γ 2 − 1 4 µ 2 . Expressing this result in terms of N (rather than µ), we obtain to leading order for large N (by replacing µ → 2N + γ, λ 2 →  For the WLβE and a general interval in the bulk, a calculation similar to that which leads to (48) yields the number variance (based on our result for β = 2 in [49], together with our conjecture (18) with [74]) andb and κb defined similarly. For the fermions in the hard box potential which can be obtained as the limit of the JβE for γ 1 = γ 2 = 1/2, we obtained the number variance for semi-infinite and finite intervals in [49] for the case β = 2. Using these results together with the conjecture (18) and its analog βπ 2 VarN In the text we have conjectured that the cumulants of order 3 and higher of the number of eigenvalues in an interval of macroscopic size in the bulk are identical for the CβE and for the GβE. This conjecture has led to predictions for fermion models. Here we provide a test of this conjecture for β = 1, 2, 4 by showing that it matches perfectly well with the rigorous results obtained for the GβE at the edge by Bothner and Buckingham [77]. The methods used in [77] and in [72,73] being completely different, this is a quite non trivial check. In Ref. [77] Bothner and Buckingham study the eigenvalues λ i of the N × N random matrices belonging to GUE,GOE and GSE ensembles near the edge, where, for large N , they scale as where the a β i form the Airy β point process. They study N [s,+∞[ the number of points a β i in the interval [s, +∞[. For v real, they prove that, for β = 1, 2, 4 and in the limit s → −∞ (i.e., towards the bulk) where k 1 = k 2 = 1 and k 4 = 2 and for β = 2 , where G (v) = log G 1 + iv π G 1 − iv π , where we recall that G(z) is the Barnes G-function [91].
which we checked explicitly using Mathematica (it is presumably equivalent to the relation (3.5) in [153]). Hence, once again, the terms of order v 3 and higher exactly coincide in the two formulae mentioned above.
Case β = 4. In that case t = v 2 √ 2π . One checks that 1/2 times the last term in (215) is log . This must be compared with the last line (214). This looks even more hopeless than for β = 1. However, using Mathematica we have discovered the identity valid for real v (where the right hand side appears to be an even function of v) whose derivation we leave as a challenge to the reader [161]. Thus, also for β = 4, the terms of order v 3 and higher exactly coincide in the two formulae mentioned above.

F Matching the variance near the edge for the harmonic oscillator
In this Appendix we show that our bulk result for the variance for the harmonic oscillator for general β, given in (40), matches for β = 1, 2, 4 the universal edge behavior obtained in [77]. The FCS formula (213) in Appendix E was obtained in [77] for the GβE. We now translate it in the context of the fermion model in the harmonic potential V (x) = 1 2 x 2 . The connection is simply a scale transformation x i = β 2 λ i , see Table 1. For general β the right edge is thus at position x + = √ βN and the width of the edge region is w N = We begin with the simplest case, β = 2. Using the expansion of the Barnes-G function [91] we find the leading terms in the expansion of (213) in powers of v: The coefficient of v 2 /2 in the expansion of (222) corresponds to the second cumulant (the variance) and therefore it gives the asymptotic behavior of the scaling function V 2 (â) from (58) forâ → −∞ as V 2 (â) 2 3 2 log (−â) + c 2 + 2 log 2 2π 2 , which, together with (58), matches exactly the bulk result (59). Similarly, for β = 1, 4 for the harmonic oscillator it was conjectured [41,95] that there exist universal scaling functions such that in the edge region Applying the correspondence (220) we obtain from (213) Again, the coefficients of v 2 /2 in these expansions give the asymptotic behaviours of V 1 (â) and V 4 (â), which, together with (223), can be summarized as which matches the bulk result (40). The leading order (logarithmic) term in (65) was conjectured in [41] for any β based on the expected matching with the bulk.