$n$-body anti-bunching in a degenerate Fermi gas of $^3$He* atoms

A key observable in investigations into quantum systems are the $n$-body correlation functions, which provide a powerful tool for experimentally determining coherence and directly probing the many-body wavefunction. While the (bosonic) correlations of photonic systems are well explored, the correlations present in matter-wave systems, particularly for fermionic atoms, are still an emerging field. In this work, we use the unique single-atom detection properties of $^3$He* atoms to perform simultaneous measurements of the $n$-body quantum correlations, up to the fifth-order, of a degenerate Fermi gas. In a direct demonstration of the Pauli exclusion principle, we observe clear anti-bunching at all orders and find good agreement with predicted correlation volumes. Our results pave the way for using correlation functions to probe some of the rich physics associated with fermionic systems, such as d-wave pairing in superconductors.

Pair correlations for photons were first considered in attempts to explain the Hanbury-Brown and Twiss (HBT) effect, where correlations between intensity fluctuations were observed for a thermal light source.[1].The increased probability of two photons arriving at a detector simultaneously compared to random chance (termed bunching) is caused by the constructive interference between individual photons.Glauber famously reconciled this effect with a full quantum description of coherence based on n-body correlation functions, which started the field of quantum optics [2].The experimental realisation of trapped neutral atoms at ultracold temperatures with large de Broglie wavelengths opened the possibility of conducting equivalent optics experiments with atoms rather than photons, termed quantum atom optics.There have since been a number of studies on the measurement of nth-order bosonic correlation functions in a variety of ultracold systems [3][4][5][6][7][8][9][10][11][12][13][14].These explorations are significant for investigating quantum statistics, quantifying information on the coherence and size of a quantum source as well as providing deeper insight into many-body quantum behaviour.
Ultracold atoms also open up the fascinating possibility of measuring fermionic correlation functions, which have no equivalent in classical optics.In the HBT effect for fermionic fields, the anti-symmetry of the wavefunction leads to destructive interference between possible propagation paths and thus a decreased probability of simultaneous detections [15].This direct demonstration of the Pauli exclusion principle is referred to as antibunching and, unlike bosonic bunching, has no classical counterpart [16].However, partly due to there being less ultarcold fermionic experiments compared to bosons [17], there have been only a handful of studies on the secondorder correlation function of neutral fermions [15,16,18] and a single measurement of 3 atom correlations [19], but nothing on higher-order fermionic correlation functions.
Here we present the first measurement of the nor-malised fermionic correlation functions simultaneously from second to fifth order for a ballistically expanding degenerate Fermi gas (DFG) of 3 He* atoms.For comparison, we also present the second-order correlation function of a thermal cloud of 4 He* atoms, which is used to sympathetically cool the 3 He* atoms, as they each display distinct behaviour [21][22][23].We utilise a Multi-Channel Plate (MCP) with a delay-line detector (DLD) to reconstruct the distribution of an ultracold 3 He* and 4 He* mixture [24] with single atom resolution in the far field after ballistic expansion [20].This allows us to directly calculate correlation functions for both bosonic ( 4 He*) and fermionic ( 3 He*) atoms, with our only limiting fac-MCP det ect or y x z g FIG. 1. Schematic showing how correlation functions are measured experimentally.An initially trapped cloud (top, blue) is released and expands as it falls onto an MCP detector (grey), which measures the arrival times of each atom as well as their x-y spatial locations.To construct the unnormalised thirdorder correlation function (G (3) (τ1, τ2)), all subsequent atoms arriving within a spatial volume ∆x and ∆y of each atom have their arrival time differences τ1 and τ2 (bottom right) recorded and histogrammed.FIG. 2. The second-order normalised correlation function g (2) (τ ) for (a) fermionic 3 He* atoms and (b) thermal bosonic 4 He* atoms.The measured correlation amplitudes for (a) and (b) are g (2) (0) = 0.84(1) and g (2) (0) = 1.11 (6), while the correlation lengths are lt = 240 (10) µs and lt = 180(40) µs, respectively.The data in (a) has ∆t = 133 µs, ∆x = 130 µm and ∆y = 560 µm, and is averaged over 2,000 experimental runs.The data in (b) uses bin widths of ∆t = 100 µs, ∆x = 130 µm and ∆y = 420 µm and is averaged over 300 experimental runs, with each run containing 116 separately out-coupled clouds of atoms, similar to [20].The errors are estimated from the standard deviation of the counts of correlated tuples across all experimental runs.
tor on maximum correlation function order achievable being data rates and detector resolution.We are able to observe fermionic antibunching for every order up to n = 5.We also investigate how the correlation width of the fermionic cloud varies with cloud size.Our results agree with theory [21,25] when considering the effects of finite detector resolution and binning size.
Correlation functions were introduced by Glauber [2] to characterise the coherence between an n-tuple of particles in space and time.We consider an nth order correlation function that considers n-fold coincidence count rates where Ψ † (r i , t i ) is the field operator for a particle at position r i at time t i and the angle brackets denote averaging.Since here we will only be studying equilibrium distributions, we can ignore the t i variable and only consider G (n) (r 1 ; . . .; r n ).To physically interpret the nth order correlation function, we take the normalised version where ρ(r i ) = G (1) (r i ; r i ) = Ψ † (r i ) Ψ(r i ) , which gives the probability of detecting n particles at points (r 1 ) to (r n ) with respect to random chance.For instance, a measurement of g (n) (0, . . ., 0) = 1 implies that the detection between the various points follows an uncorrelated Poissonian distribution.In contrast, g (n) (0, . . ., 0) > 1 implies there is bunching present in the system and g (n) (0, . . ., 0) < 1 implies anti-bunching.
Our experiment starts with a combined degenerate Fermi gas of 3 He* atoms and a Bose-Einstein condensate (BEC) of 4 He* atoms, both trapped in a magnetic trap with trapping frequencies for the 3 He* atoms of ω x,y,z = 2π × (58(3), 694(1), 701(2))Hz [24,26].After the trap switches off, the clouds fall ∼850mm (fall time t T OF = 416ms) onto a MCP and DLD, which can detect the 3D location of individual atoms with ∼130µm x, y and ∼3µs z resolution [20,27].Due to the different masses of the two helium species, applying a small magnetic field gradient during time-of-flight (TOF) separates the arrival times of the clouds, allowing the distribution of each species to be measured in a single shot [24].The TOF is sufficiently long that the detected position distribution at the detector is close to the cloud's in-trap momentum distribution and hence the correlation functions correspond approximately to the in-trap momentum.
The single atom resolution in 3D allows us to reconstruct correlation functions directly.For simplicity, we only consider correlations along the single axis with the highest detector resolution (the z axis of the ballistic expansion).Since this corresponds to the arrival time of a falling cloud at the detector, we refer to this as correlations in arrival time t.We can therefore consider the simplified volume integrated correlation function, which for the second order is Here, we have averaged over all arrival times and are hence now finding the correlation between detected events with a τ delay between their arrival times.The area of the spatial integral over the x and y directions is chosen to be comparable to the correlation lengths in these directions so that all events within the integration volume are correlated.We will also extend this logic to the nth order correlation to obtain the simplified volume integrated correlation function g (n) (τ 1 , . . ., τ n−1 ), as shown in Fig. 1 for g (3) (τ 1 , τ 2 ) (see [25] for further details).In Fig. 2 (a) and (b), we show the experimentally measured g (2) (τ ) for a cloud of ultracold bosons (a) and FIG. 3. The experimentally measured temporal correlation length lt (a) and correlation amplitude g (2) (0) (b) for various time-of-flight widths of the 3 He* cloud.These are compared to the theoretically ideal value (dashed line) and expected value taking into account experimental factors (solid line) where realistic effects such as finite binning, resolution and dark counts are included via a Monte Carlo simulation [25].For (b), the theoretically ideal value is g (2) (0) = 0 for all cloud widths, and hence it is not included in the plot.We find excellent agreement for both correlation length and amplitude between the experimentally measured values and the expected non-ideal values.
By fitting Gaussians to the data in Fig. 2 (a) and (b) [25], we find the correlation length at the detector for the DFG is l t = 240(10) µs, while for the thermal 4 He* cloud it is l t = 180(40) µs.Incorporating experimental factors such as finite binning and detector resolution (as in previous investigations [6]), we perform a Monte Carlo simulation (see [25] for details) and find an expected correlation amplitude of g (2) (0) = 0.83 for the DFG, which agrees well with the experimentally measured value of 0.84 (1).We similarly find the expected correlation length for the 3 He* atoms l t = 269 µs, which again agrees well with the measured value of 270(10) µs.
To further investigate the dependence of the antibunching on temperature, we vary the final temperature of the 3 He* atoms by evaporating more of the coolant gas (bosonic 4 He*).This reduces the temperature while keeping the total number of fermions approximately the same.As there are few bosons left in the mixture and most of those remaining are in the condensate, it becomes difficult to ascertain an accurate measure of temperature from the 3 He* cloud, as discussed in Ref. [24].Hence, we use the time-of-flight width of the Fermi cloud, given by σ T OF = t T OF γ(ξ) k B T m (see [25] for details), as a proxy for temperature [28].Note that the time-of-flight width is isotropic even for a degenerate Fermi cloud.Fig. 3 shows the measured correlation length l t (a) and maximum anti-bunching amplitude g (2) (0) (b) as a function of σ T OF , along with theoretical expected values [25].
In conclusion, we have measured the fermionic normalised correlation functions from second to fifth order, with orders greater than three being presented for the first time.The observed anti-bunching effect in n-tuples of 3 He* atoms agrees with the predictions of Wick's theorem combined with realistic experimental effects.Our method could be extended to in principle measure arbitrarily high-order correlations, with the major limiting factor being data acquisition rates.The ability to measure higher-order correlation functions of fermionic atoms opens up a range of fascinating possibilities with atom number estimate.See [24] for further details.
[29] Correlation function are normalised in such a way as to attain unit value in the case of perfect coherence.For the first order correlation function this means the normalisation differs slightly and is given by g (1) (r1, t1; r2, t2) = , with the volume integrated version given by g (1) (τ ) = dr dt G (1) (r,t;r,t+τ ) .
[30] To clarify this point further it is easier to consider the form g (n) (0) = ν∈P (n) η i(ν) g (1) (0) |supp(ν)| (see [25] for derivation) where P (n) is the set of permutations of the numbers 1 to n, supp(ν) is the support of the permutation ν, and i(ν) is the number of inversions contained in the permutation.Given that there is an equal number of permutations with odd and even numbers of inversions, we will have an even number of positive negative, which will exactly cancel to give zero.[35] The support of a permutation is the set of numbers (or more generally elements) that are moved from their original positions under that permutation.For example, the support of the permutation 1234 → 2134 is {1, 2}.

SUPPLEMENTARY MATERIAL
Decomposition of g (n) via Wick's theorem The general form of Wick's theorem can be simplified for a correlation function using its respective anti-commutation or commutation relations (see example 5 in [34]) Here P ([1, n]) is the set of permutations of the numbers 1 to n, i(ν) is the number of inversions contained in the permutation ν, Nk = Ψ † k Ψk is number operator, ν(j) represents the jth term in the permutation ν and η is +1 for bosons and −1 for fermions.We can readily identify the LHS of Eqn. 6 as G (n) (r 1 , t 1 ; . . .; r n , t n ) (for now we will consider G (n) (t 1 ; . . .; t n ) for notational convenience and relevance, however, this derivation is also valid for an arbitrary position vector in some parameter space) and the RHS as a function of various G (1) 's.Using the notation of time differences, we can then rewrite Eqn.6 as where G (1) (t k ; t m ) = Ψ † j Ψν(j) by definition.As we wish to consider the normalised correlation functions, we divide both sides by ρ 1 . . .ρ n to give Note that we have used the fact that each term in the sum has exactly one Ψ † k and Ψk for each k between 1 to n and have 'assigned' each of these a normalisation of √ ρ k from ρ 1 . . .ρ n .If ν(j) = j we have G (1) (t j ; t j ) = ρ j , hence by definition, these terms are cancelled by the denominator, otherwise we obtain g (1) (t j ; t ν(j) ).Thus, g (1) (t j , t ν(j) ), (9) where supp(ν) is the support of the permutation ν [35].We can use Eqn. 9 to derive any order correlation function.
Tab.I shows the correlation functions up to n = 5.
To match the measured correlation function, we simplify further by substituting t k = t + τ k−1 and integrating over t to obtain the averaged correlation function, dt g (1) (t + τ j ; t + τ ν(j) ) (10) where we have shifted the indexation down by 1 to reflect the relabeling and τ 0 = 0. Next notice dt g (1) (t + τ j ; t + τ ν(j) ) = dt ′ g (1) (t ′ ; t ′ + τ ν(j) − τ j ), thus we can write these terms as g (1) (τ ν(j) − τ j ).From this, we see the exact form of the decomposition of the average nth-order correlation function, g (n) (τ 1 , . . ., τ n−1 ), in terms of g (1) (τ )'s, is given by Wick's theorem as g (1) (τ ν(j) − τ j ).(11) In the text we consider the diagonal correlation function g (n) (τ ) ≡ g (n) (τ, ..., τ ), and specifically we want to understand the behaviour of g (n) (τ → 0).Note τ can be positive or negative, however, we assume the same sign for all particles (i.e.τ 1 , . . ., τ n−1 all have the same sign).In the ideal case of g (1) (τ → 0) = 1, Eqn. 11 gives the interesting result that if n ≥ 2 and η = −1 then g (n) (τ ) = 0 for all τ .This intuitively makes sense, as g (n) (τ ) > 0 for n ≥ 2 would require two fermions to be measured at the same location, which in the ideal case is impossible.We will TABLE I. Decomposition of the nth order fermionic (η = −1) correlation functions, for n from 2 to 5, in terms of the first-order correlation function.The notation S k represents the set of k element subsets of the set S. Alternatively, this can be thought of as the set of possible outcomes for the operation S choose k.We have also used the notation DR(S) = {k ∈ M (D(S)) : |k| = R}, i.e. the set of elements of (D(S)) with cardinality equal to R, where D are the derangements of the set S b and M acts on a set of permutations X as follows M (X) = {i, ν(i)} : i ∈ [1, |ν|] : ν ∈ X , giving the mapping set b .Generally, ti can be thought of as some arbitrary dimension position vector.For comparison, we show the simplified decomposition of g (n) (0), where all ti = 0.

2R
(i,j)∈k g (1) (ti, tj) a Permutations of S such that no element remains at the same place.This can also be stated as the support of the permutation is equal to the original set S. b This can be thought of as the set of the sets of pairs of elements that map to each other for each permutation contained in set X. d In this column of the table, we have used the indexation shorthand (i, j) ∈ S for {{i, j} : {i, j} ∈ S and i < j}.The purpose of the implicit condition i < j is to disambiguate the indexation as if {i, j} ∈ S then {j, i} ∈ S. Thus if a sum or product over {i, j} ∈ S were followed strictly it could lead to double counting.
see, however, there is deviation from this ideal case experimentally due to various effects.Hence, to predict the value of g (n) (0) under realistic measurement conditions we set τ = 0 (assuming nothing of the value of g (n) (0)) and reach This can be rewritten as where T (k) =!k for bosons (η = 1), with !k denoting the sub-factorial of k, and T (k) = (k − 1) for fermions (η = −1).
To understand this simplification consider that for |supp(ν)| = k, we must displace k elements of our original set of n, so we n k choices of the set of elements that we are displacing.For bosons we then simply consider the number of possible ways we can permute these k elements such that none are left in their original positions, known as a derangement, of which there are !kways to do for such a number of elements.For fermions, we instead need to consider how many derangements have a different parity of inversions.This is given by (−1) k−1 (k − 1), noting this difference is independent of the original set of elements we chose.
For g (1) (0) = 1 we obtain the expected relations of g (n) (0) = n! for bosons and g (n) (0) = 0 for fermions.For bosons, this is can be understood by noting that each term in the sum from Eqn. 12 contributes 1, and thus the total sum equals the number of terms, i.e. the number of permutations of the numbers from 0 to n − 1, which is n!.For fermions, we see that the sum in Eqn. 12 is equal to the difference between the number of permutations with odd and even numbers of inversions, which is always equal, and hence we obtain 0. As discussed in the main text, this implies that the measured value of g (n) (τ → 0) can be predicted solely by the experimentally measured value g (1) (0).It can be seen that the measured value rapidly deviates from zero for g (1) (0) ̸ = 0, explaining the measured g (n) (τ ) distributions (see main text Fig. 4).
Exact form of correlation length for non-interacting fermi gases in a harmonic trap In the text while we are measuring momentum correlations, we construct these correlations using a time-of-flight distribution and hence discuss correlators in terms of detected position and time, however for this derivation we will use the equilibrium position and momentum of the trapped atoms, i.e. the underlying distribution which produces the time-of-flight profile.We discuss how to convert these momentum correlation lengths to their time-of-flight equivalents in the text Sec. .From Sec. , we can see that if we find an analytical form for the correlation length of G (1) (p 1 ; p 2 ) we can use Wick's theorem to extend it to find an analytical form for any G (n) (p 1 ; . . .; p n ).The definition of the first-order coherence function can be rewritten in terms of the Wigner function W (p, q) [21], G (1) with r, r ′ and p, p ′ referring to the in-trap position and momentum, respectively.Under the local density (or semiclassical) approximation, which assumes W (p, q) is equivalent to a spatially homogeneous system with a varying chemical potential, the Wigner function of harmonically trapped atoms is given by where α = m 2k B T , β = 1 2mk B T , q = (ω x x) 2 + (ω y y) 2 + (ω z z) 2 is the normalised position vector, and ξ = e µ/k B T is the fugacity.Notice that exchanging the momentum p and position q in a harmonic trap is equivalent to exchanging the values of β and α.Thus, the momentum and position correlation functions are functionally equivalent for noninteracting gases in a harmonic trap.
The main quantity of interest measured in this work is the average normalised correlation function.We consider the full definition g (1) (∆p) ≡ dP G (1) (P − ∆p 2 ; P + ∆p 2 ) dP ρ(P − ∆p 2 ) ρ(P where ρ(p i ) = G (1) (p i ; p i ) is the momentum density function.Using Eqns.15 and 16, we can approximate the solution of Eqn. 17 (assuming ∆p ≪ (2mk B T ) 2 ) as with correlation length l i = ℏ 2si γ(ξ), where is the size of a thermal cloud in the i-axis of a harmonic trap and γ(ξ) = Li3(ηξ) Li4(ηξ) .This approximation accurately predicts the bulk behavior of the correlation functions, even in highly degenerate regimes.Accordingly, we can express g (2) (∆p) and g (3) (∆p 1 , ∆p 2 ) for fermions as Gaussian functions with corresponding correlation lengths over the whole temperature range of interest.Following from the relations in Tab.I, we find Note that the sign of the ∆p 2 + ∆p 1 terms is due to us considering an ordered set of events.If we were looking at an unordered set, we must consider an average of ∆p 2 + ∆p 1 and ±(∆p 2 − ∆p 1 ).

Predicted experimental correlation amplitudes and lengths
While the results from Sec. and indicate the expected correlation length and amplitude at a given order for a perfect detection system, experimental measurements will typically yield reduced amplitudes, due to a variety of effects.To understand the most relevant effects, we will briefly describe how the nth order correlation function is computed for this particular work.We first compute the spatial (both in x and y axes) and temporal differences between all pairs of atom detections of a given atomic species.For the nth order correlation function, we consider all n-tuples of particles, for example, pairs in second order and triplets in third order, whose spatial difference between all pairs of particles in the tuple are less than a given threshold ∆x and ∆y for the x and y axes respectively.This is practically equivalent to integrating the full correlation function over these dimensions.Hence, we expect improved correlation amplitudes for decreasing ∆x and ∆y at the cost of signal-to-noise performance.We then take all the temporal difference τ 1 , . . ., τ n−1 , repeating for all counts in a given experimental realisation.These differences can then be histogrammed to obtain the unnormalised n-th order correlation function along the time dimension.The unormalised distribution is averaged over many experimental realisations to improve signal-to-noise.To normalise the distribution, we repeat the above procedure, however this time for counts recorded in different experimental realisations, which cannot possibly interfere and hence should contain no underlying correlation allowing us to determine the bulk distribution.Note that detection efficiency does not affect the normalised correlation function, as it proportionally affects the numerator and the denominator by the same amount.From this, we see that the strongest deviation from the ideal case comes from the effective integration over the x and y axes due to the finite bin size.There is also a slight integration over the time axis due to again finite bin size.This effect compounds at higher orders, as for each new particle added we integrate over another set of axes.Less noticeable, but still significant effects are the resolution of the detector and the dark, or background, count rate of the detector.
To account for these effects, we can perform a Monte Carlo simulation of our full measurement process, as described above, assuming the underlying true distribution follows the relations in Sec. and and with the various imperfections included.As our detection parameters, such as resolution and dark count rate, are fixed, and the ideal correlation amplitude is always zero for any order, the only relevant variable in our simulation is the correlation length.The correlation length of the trapped gas, as derived in Sec., is given by l i = ℏ 2si γ(ξ).As mentioned above, in practice, we experimentally measure the correlations in a ballistically expanding cloud.Hence, we rescale the momentum space correlation length to give the experimentally measured TOF correlation lengths, where t T OF is the flight time, and now i ∈ (t, x, y).For the t-axis, as we measure arrival time rather than a location, the TOF scaling is slightly different.If the expansion of the cloud is negligible compared to the fall distance, we can consider a correlation time of where g is the acceleration due to gravity.
We can now test the validity of our model by varying the TOF correlation length and comparing the results to the measured amplitude and correlation length.In practice, we vary the correlation length by changing the The predicted correlation length and amplitude for g (2) (τ ) are l SIM t = 250(1) µs and g (2) (0)SIM = 0.853(1) respectively.These are extracted by fitting a gaussian to the simulation output (blue line).For comparison, we also show the ideal distributions, plotted using Eqns.18 and 19, for the first and second orders, respectively (red dashed line).
temperature of the trapped cloud.However, as it is difficult to ascertain an accurate measure of temperature from the TOF profile of the 3 He* cloud alone, we instead use the (isometric) time-of-flight width of the cloud, given by σ T OF = t T OF γ(ξ) k B T m , as a proxy for temperature.From Eqn.21, we can see that l T OF i = ℏωi 2mσ T OF t 2 T OF , and thus, as ω i is fixed σ T OF alone is a sufficient initial condition for our Monte Carlo simulation and can be directly compared to experiment.
The relevant experimental parameters used in our simulation are listed in Tab.II.The resultant g (1) and g (2)  distributions for correlation lengths of (l t , l x , l y )=(244.6 µs, 83.0 µm, 996.6 µm) (σ T OF = 8 mm) with and without the experimental imperfections are shown in Fig. 6.The predicted values for the experimentally measured amplitude and correlation length for various cloud sizes are shown in Fig. 7, showing good agreement between theory and experiment and demonstrating that it is crucial to include these realistic imperfections in order to obtain accurate predictions.

3
He* atoms.Such experiments include probing interesting many-body phenomena such as time crystals [31] and d-wave superconductors [32] or investigating foundational quantum mechanical effects such as the weak equivalence principle [33].

FIG. 7 .
FIG. 7. The experimentally measured temporal correlation length lt (a) and correlation amplitude g (2) (0) (b) for various time-of-flight widths of the 3 He* cloud compared to the theoretically ideal value (red dashed line) and the results of our Monte Carlo simulation (black solid line).
Bin size (∆t, ∆x, ∆y) 100 µs, 130 µm, 420 µm TABLE II.Detection and analysis parameters used in Monte Carlo simulation.Analysis parameters are the same as those used to analyze the experimental data.