Natural orbitals and their occupation numbers for free anyons in the magnetic gauge

We investigate the properties of natural orbitals and their occupation numbers of the ground state of two non-interacting anyons characterised by the fractional statistics parameter $\alpha$ and confined in a harmonic trap. We work in the boson magnetic gauge where the anyons are modelled as composite bosons with magnetic flux quanta attached to their positions. We derive an asymptotic form of the weakly occupied natural orbitals, and show that their corresponding (ordered descendingly) occupation numbers decay according to the power law $n^{-(4+2\alpha)}$, where $n$ is the index of the natural orbital. We find remarkable numerical agreement of the theory with the natural orbitals and their occupation numbers computed from the spectral decomposition of the system's wavefunction. We explain that the same results apply to the fermion magnetic gauge.


I. INTRODUCTION
Quantum statistics describes how the wavefunction of a multi-particle quantum system changes under particle exchange.In three dimensions, particles can have statistics that are either fermionic or bosonic.Anyons are class of particles that exist only in two dimensions with fractional statistics, between that of bosons and fermions, allowing their states to gain arbitrary complex phases under exchange [1][2][3][4][5].Such a statistics plays a role in understanding the fractional quantum Hall effect (FQHE), with emergent quasiparticles that have been identified as anyons [6].In this model, anyons are assumed to form a non-interacting gas (reviewed in Section II of this paper) whose free Hamiltonian is unitarily equivalent to a system composed of either fermions or hard-core bosons with magnetic flux quanta bound to their positions [7][8][9].This picture is usually refereed to as the magnetic gauge.The magnetic flux quanta realise the exchange phases as the Aharonov-Bohm phases.Anyons also find application in quantum computing, as the exchange of non-abelian anyons (described by multicomponent wavefunctions) allows the implementation of topologically protected quantum gates which are intrinsically robust against local noise [10][11][12][13].Such anyonic quasiparticles can be detected experimentally [14] using interferometry [15][16][17], correlation measurements in collider geometries [18] or spectroscopy [19].
Due to the presence of anyonic quasiparticles in the FQHE, there has recently been interest in developing the Kohn-Sham density functional theory (KS-DFT, see [20,21] for more background) for the treatment of confined non-interacting many-anyon systems.In particular, KS-DFT has been developed for anyons described as composite flux-fermions [22,23].The reduced density matrix functional theory (RDMFT [24]) is suitable for * tomasz.maciazek@bristol.ac.uk describing strongly correlated systems [25].As anyonic quasiparticles appear in strongly correlated systems, this method has been applied to an FQHE system [26] (specifically, to a confined quantum dot in the spin-frozen strong magnetic field regime).Recent developments in the foundations of RDMFT for bosons [27] suggest that RDMFT might be a suitable tool also for studying anyon gas in the boson magnetic gauge.
In many methods of quantum chemistry the choice of the appropriate single-particle basis is a crucial step.For instance, the Hartree-Fock method approximates the Nbody wavefunction of a system as a single Slater determinant.This approach has been applied to study the fermion-flux composites in relation to superconductivity and the presence of the energy gap in the free anyon gas [7,28].More generally, the multi-configurational self-consistent field method (MCSCF) relies on superposing several Slater determinants coming from a given single-particle basis.For a given N -particle (bosonic or fermionic) quantum state, there exists a distinguished multi-configurational expansion called the natural expansion.The single-particle basis that is used in the natural expansion consists of the natural orbitals (NOs) that are defined as eigenfunctions of the state's one particle reduced density matrix (1RDM).It is known that the natural expansion has the fastest convergence of all the possible expansions [29].The eigenvalues of the 1RDM are called the natural occupation numbers (NONs) and denoted by ν n , n = 1, 2, . . ., where ν 1 ≥ ν 2 ≥ . . . .In this paper, we employ the convention for the 1RDM to be always of trace one.Due to normalisation, the NONs (ordered descendingly) have to decay to zero if the singleparticle basis is of infinite dimension.The rate of their decay reflects certain fundamental features of the system at hand.For instance, only the first N NONs of an uncorrelated state of N electrons (a Slater determinant) are nonzero.On the other hand, a slow decay of NONs characterises a highly correlated quantum state.Furthermore, the knowledge of the NOs and NONs allows one to compute the expectation value of any one-particle observable where the individual contributions of the NOs are proportional to their corresponding NONs.For these fundamental reasons, the asymptotic rate of decay of NONs in different quantum systems has been an object of notable interest in quantum chemistry.
One of the early results in this area is Hill's asymptotic [30] which concerns ground states of two-electron systems confined in an external potential with spherical symmetry.Due to this symmetry, the 1RDM has a blockdiagonal structure where the blocks are labelled by the angular momentum quantum number l. Hill's asymptotic states that the total occupancies of the blocks, ω l = n ν nl , satisfy where C H is a constant that can be computed explicitly from the given quantum state.This asymptotic behaviour of ω l was anticipated by the preceding numerical calculations for the helium atom [31,32].Numerical calculations for the harmonium atom [33,34] further confirmed the validity of Hill's asymptotic beyond Coulomb external potentials.The large-n asymptotic of ν nl (for fixed l-sector) has proved to be a more difficult problem to study due to the lack of sufficiently accurate electronic structure data.Various conjectural forms of this asymptotic have been proposed and subsequently refuted over the years [32,35,36].Finally, its correct form has been determined for ground states of two-electron systems in central potentials and proved rigorously to be [37] where the constant A does not depend on l and can be calculated explicitly from the wavefunction at hand.This result has been subsequently generalised to singlet states of systems with arbitrary symmetries and numbers of electrons [38].Another work concerning general quantum systems of N Coulomb-interacting particles proved the asymptotic [39] lim where the constant B can also be calculated explicitly from the wavefunction [38,39].Note the lack of the subscript l in the formula (3) that does not assume any symmetries.The above results have subsequently led to the recent discovery of a universal power law governing the accuracy of wave function-based electronic structure calculations [40].Similar results have been recently obtained for a system with the Fermi-Huang interparticle interaction [41].
For the reliable application of quantum chemistry methods to anyonic systems it is necessary to understand similar asymptotic behaviour for models of noninteracting anyons.Here it is relevant that anyonic systems are two-dimensional as opposed to the previously mentioned quantum-chemical systems which are threedimensional.Due to this change in dimensionality, new technical tools have to be applied in order to extend the methodology of the above cited papers to anyonic systems.However, the core feature remains true: the NOs and NONs are solutions to the eigenproblem of an integral operator whose kernel has a particle-particle coalescence cusp that drives the large-n asymptotic of the NOs and NONs.More specifically, as we explain in Section II, the coalescence cusp is proportional to |z 1 − z 2 | α , where z 1 , z 2 are the positions of the two boson-flux composite particles (represented by complex numbers) and α ∈ [0, 1] is the fractional statistics parameter.In consequence, we show that the natural amplitudes (NAs, the positive square roots of the NONs, σ 2 nl = ν nl ) of the ground state of the system at hand satisfy and provide explicit expression for the constant D(α) in Equation (37) in Section III.We also find very accurate asymptotic forms of the weakly occupied NOs in terms of appropriately transformed integer Bessel functions.

II. THEORETICAL BACKGROUND
Anyons have quantum statistics between fermions and bosons -the exchange of a pair of abelian anyons results with the multiplication of the many-anyon wavefunction by the phase factor e iπα , where α ∈ [0, 1].Consequently the wavefunction Ψ α of N abelian anyons can be expressed in terms of a bosonic wavefunction Ψ B as where the complex number z j = x j + iy j describes the position of the j-th particle, and α is the fractional statistics parameter (see e.g.[7,42]).Similar mapping defines anyons in one dimension, where the one-body reduced density matrix of anyons confined in a harmonic trap and its NOs can be computed efficiently even for large particle numbers [43], showing that the largest NON follows a power law ν 1 ∼ N p(α) with 0 < p(α) < 1.The anyonic wavefunction Ψ α interpolates between bosons for α = 0 and fermions for α = 1.The free-particle Hamiltonian Ĥfree = ∇ 2 1 + ∇ 2 2 acting on Ψ α then transforms under the gauge transformation in Equation ( 5) to a magnetic Hamiltonian acting on Ψ B [2,28].In the case of two particles, this magnetic Hamiltonian in centre of mass coordinates can be expressed as [44] where Z is the centre of mass, z = z 1 − z 2 is the relative coordinate and A(z) is the magnetic vector potential where eΦ = αℏ.In this model, the anyonic exchange phases can be thought of as Aharonov-Bohm phases due to the presence of a magnetic vector potential in the relative Hamiltonian.The bosonic Hamiltonian therefore describes boson-flux composites with hard cores.The eigenstates for this two-anyon system in a harmonic potential are known [44].Under the potential V external = (z 2 1 + z 2 2 )/2 (here and in the following the "atomic" units are used in which the harmonic potential strength, the particle mass and ℏ are all equal to one) the ground state wavefunction in the boson-flux composite picture is given by The natural occupation numbers (NONs) ν n for a many-particle quantum system are defined as the eigenvalues of the one-particle reduced density matrix (1RDM).Recall that for a two-particle wavefunction in R d the 1RDM reads [45] and the natural orbitals ϕ n are its eigenstates, i.e.
If the wavefunction Ψ is real and symmetric, it is known that its NOs and NONs can be also found by solving the following homogeneous Fredholm equation of the second kind [46], which circumvents the necessity of computing γ In the above equation the eigenvalues {σ n } ∞ n=1 are the natural amplitudes (NAs), ν n = σ 2 n .Moreover, if the wavefunction Ψ is rotationally symmetric (as is the case for any eigenfunction of the Hamiltonian ( 6)), then the relevant eigenproblem is blockdiagonal where the blocks are enumerated by the angular momentum quantum number l.In summary, the rest of this paper will be devoted to asymptotically solving the integral equation where B is given by formula (8), for an arbitrary fixed spin sector l, in the limit of n → ∞ which means studying the weakly occupied NOs and their corresponding NAs.The key intuitions here are that i) the radial part of ϕ nl is highly oscillatory (as shown in Fig. 4a-c -the nth numerically computed NO has n − 1 nodes), ii) when integrating a highly oscillatory function against a function that has discontinuous derivatives, the result goes to zero as a polynomial of the inverse of the oscillation frequency.

III. DERIVATION OF ASYMPTOTICS
Restricting to a fixed-l sector means taking the NOs of the following forms where the factor e −r 2 /2 is included for the sake of convenience.The orthonormality condition of the NOs then reads Using the relative angle variable θ 12 = θ 1 − θ 2 and the polar coordinates z j = r j e iθj , j = 1, 2, the Fredholm Equation ( 12) becomes In order to evaluate the angular integral above, note that the function cos (θ 12 l) can be written as a polynomial of degree 2l in the variable t where T l is the Chebyshev polynomial of order l, which then allows the finite polynomial expansion where the coefficients a (l) p can be found using the explicit expansion formulas for T l .In particular, Plugging this expansion into the LHS of ( 15), we get Next, we use the result that for any β ∈ R, 2π 0 where and 2 F 1 is the ordinary hypergeometric function.
So far, all the calculations have been exact.However, we will next start making approximations in order to make the large-n asymptotic solution to the Equation (15) more tractable.The key fact to notice is that the leading contribution to the n → ∞ asymptotic of σ n comes from the lowest order cusp of the LHS of Equation (15) around r 1 = r 2 .Let us next look at the asymptotic expansion of G β (r 1 , r 2 ) around r 1 = r 2 when 0 < β < 1.The only contribution to the cusp around r 1 = r 2 is the first term in the RHS of Equation (19).Its asymptotic expansion reads (20) The first three terms of this expansion are plotted in Fig. 1.In the above formula, we can see that the leadingorder cusp is |r 1 − r 2 | β+1 .Similarly, when β > 1 the analogous cusp will appear in the order β + 1.This in turn means that the leading-order cusp in the Equation (17) will come from the (p = 0)-summand.Thus, the large-n asymptotic solution to the Equation ( 15) can be obtained from the following equation which extracts only the leading-order cusp around (21) In the above equation we have also used the fact that a (l) 0 (r 1 , r 1 ) = 1.
In analogy to the three-dimensional case [37] we choose an ansatz for ψ n,l of the form where Functions J l are the Bessel functions of the first kind.
As previously mentioned, taking the limit n → ∞ means considering highly oscillatory NOs, thus we necessarily have κ nl >> 1.The strategy is now to extract the leading order expansion of the Equation ( 21) in the powers of 1/κ nl which will lead us to certain consistency condition for the function g nl (r).Because κ nl is large, we can approximate the Bessel functions as [47] J For convenience, we will ignore the phase shift proportional to π/4 under the cosine, as it will not alter the resulting consistency relations.Consequently, the integral equation (21) becomes g nl (r 1 ) cos (κ nl g nl (r 1 )) .
With the change of variables Next, we extract the leading asymptotic around u 1 = u 2 using the fact that This allows us to write the integral Equation (24) as where u ∞ = lim r→∞ g nl (r).Finally, we use a theorem from the paper [48] to extract the leading order behaviour when in κ nl → ∞ (see Appendix A for more details).This leads to Equating coefficients of cos (κ nl u 1 ) and expressing everything back in terms of the variable r, we arrive at the following consistency condition defining g nl (r) where erf is the error function, and In order to determine the dependence of σ nl and κ nl on n together with the form of the function f nl , we refer to the orthonormality condition of the NOs ( 14) which takes the explicit form In order to satisfy the above equality, we will make use of the orthogonality relation for Bessel functions [47] where µ n is the nth node of J l .Note first that if we require the relation to be satisfied and requiring that g n1l (r) = g n2l (r) ≡ g l (r) we can readily apply the identity (31) to the LHS of (30) after the familiar change of variables u = g l (r) under the integral.Additionally, recall that the nth NO must have n − 1 nodes.In light of Equation ( 27), this requires identifying the functions g nl (r) and κ nl as Consequently, we get that From the relation (32) we determine that f nl (r) takes the form where the normalization factor reads Using the well known asymptotic for the nodes of the Bessel functions µ n ≈ nπ, we summarize the results of this section as follows.

IV. NUMERICAL METHODS
In this section, we carry out the numerical verification of the asymototics predicted in Equations ( 37) and (38).This is done by writing the Equation (12) in the (κ-scaled) harmonic oscillator eigenbasis where κ > 0 and L |l| m is a generalised Laguerre polynomial.These calculations are specific to l = 0.As will be explained later, the parameter κ will be optimised which allows us to pick the the basis in which the NONs converge at the fastest rate.To simplify the notation we omit the l-subscripts, i.e. we write φ(κ) m,0 ≡ φ(κ) m .In the truncated basis ϕ (κ) m , m = 0, 1, . . ., M the Equation (12) for l = 0 is solved via the diagonalisation of the (M + 1) × (M + 1) matrix Let us next briefly describe the steps that we apply in order to evaluate some of the multidimensional integrals in the LHS of Equation (40).Similarly to the methodology of Section III, we use the relative angle θ 12 and compute the corresponding integral over θ 12 using the result from Equation (18).Next, we change the the radial coordinates r 1 , r 2 to r 1 = r cos ξ and r 2 = r sin ξ, 0 ≤ r ≤ ∞, ξ ∈ [0, π/2].The Jacobian of this transformation is r.
Under this change of variables, we have where the function K α (ξ) reads (after using an identity for the hypergeometric function to simplify its form) Next, we apply the polynomial expansions of the Laguerre polynomials.It turns out that the integrals over r can be computed analytically in terms of the Euler gamma functions.Thus, we are only left with the task of evaluating numerically the integrals over ξ.After the above transformations, the matrix elements read where the integration over ξ appears only in the integrals Note that in order to evaluate the expressions (42) numerically for all 0 ≤ m 1 , m 2 ≤ M , we only need the integrals J (α) k for k = 0, 1, . . ., 2M which can be precalculated separately.The difficulty of this approach is that we need to know the values of the integrals J (α) k with very high precision, because they enter the alternating sums in Equation (42).To this end, we have used P ython's library mpmath.For the calculations presented in this section, we have set the size of the one-particle basis to M = 400 and the precision to 2M .Note also that the expressions (42) can be computed efficiently for all 0 ≤ m 1 , m 2 ≤ M at once using the vectorisation technique, i.e. by recognising that the Equation ( 42) allows one to express matrix A (α) (κ) as a product of three matrices.
The optimal paramater κ has been chosen by maximising the fidelity It has turned out that the optimal value of κ is approximately independent of α and for M = 400 it is κ opt ≈ 7.
The resulting NAs are plotted in Fig. 2. Thanks to the optimal choice of the parameter κ, we have obtained the convergence of the first 105 -125 highest NAs (the exact number depends on α) with the single-particle basis of the size M = 400.In Fig. 2 we have plotted the values of n α+2 σ n vs. 1/n, which allowed us to determine the values of the constants D(α) with the relative accuracy of the order 10 −3 .As shown in Table I, this is in perfect agreement with the theoretical values calculated from Equation (37).
In Fig. 4 we have plotted the numerically calculated NOs and compared them with their asymptotic forms from Equation (38).We observe remarkable agreement of the asymptotic form even for values of n as low as 30.For the convenience of comparison, in the bottom row of Fig. 4 we have plotted the values of ϕ n (r) √ re r 2 /(2(α+2)) vs. g(r) which extracts the oscillatory part of the NOs.One can see that the NOs converge to their respective asymptotic forms very fast.

V. CONSEQUENCES OF THE NO AND NA ASYMPTOTICS FOR COMPUTING ANYON CORRELATIONS
The asymptotic properties of NOs and NAs derived in the preceding sections are of broad relevance to physics of anyonic systems astthey show universal properties of anyonic wavefunctions.In particular, our results show that when an anyonic wavefunction is expanded in a single-particle basis, the convergence of such an expansion is much slower than in the case of fermionic or bosonic systems.This has fundamental consequences for applications of any numerical methods to anyonic systems.These consequences are particularly evident in two-particle systems where the natural orbitals and natural occupation numbers allow one to reconstruct the two-anyon wavefunction (up to a diagonal rotation of the single-particle basis), making it possible to compute the expectation value of any quantum observable.Such twoanyon systems have been extensively studied from the point of view of anyon correlations and anyon interferometry [49][50][51][52][53][54] (including recent experimental proofs of anyonic statistics [16,18]), where one of the relevant observables is the two-body observable called the bunching parameter [54].
In this section, we show in detail how to calculate such expectation values using the example of two-particle correlation coefficient τ defined as where we parametrise the configuration space of the anyons via z j = r j e iθj , j = 1, 2 with θ 12 = θ 1 − θ 2 being the angle between the position vectors of the anyons.The correlation coefficient τ has been originally introduced in the context of electronic wavefunctions [55], however it is closely related to the anyon bunching parameter introduced more recently [54].
In order to calculate the correlation coefficient τ for the two-anyon system at hand, we will work in the boson magnetic gauge, i.e. use the wavefunction Ψ (α) B defined in (8).It is convenient to change the variables to the centre of mass, Z, and relative position, z, given by In the new coordinates, the two-anyon wavefunction has the product form What is more, the correlation coefficient also takes a simple form in the new coordinates, i.e.
It is a matter of a straightforward calculation to find the relevant expectation values.The result reads Recall that for α = 0 the quantum state B is just a product state of two bosons which is uncorrelated.Consequently, the correlation coefficient vanishes for α = 0. What is more, τ is a monotonic function of the statistics parameter α and it reaches its minimum value τ = −1/3 for α = 1, i.e. when the wavefunction Ψ (α) B describes a strongly correlated state of two hard core bosons.
Let us next take a closer look at how to compute the above correlation coefficient τ with the NOs of Ψ (α) B .The two-anyon wavefunction can be written as Note that only those terms for which |l 1 − l 2 | = 1 contribute to the numerator; otherwise ⟨ϕ n1,l1 | z |ϕ n2,l2 ⟩ vanishes.Although the closed forms of the expectation values in the formula (51) are difficult to find, numerics using the approximate NOs (38) shows that the convergence of the expression ( 51) is rather slow, requiring the calculation of thousands of terms to get numerically close to the analytic value (49).

VI. DISCUSSION AND CONCLUSIONS
In this work we have analysed the asymptotic behaviour of the natural orbitals and their occupation numbers/natural amplitudes in the ground state of two noninteracting anyons in the boson magnetic gauge.We have derived exact asymptotic forms of the natural orbitals from a fixed l-sector of the 1RDM when the ordinal number of the orbital n that indexes NOs and NAs tends to infinity (the latter are arranged descendingly according to their absolute values).While the asymptotic forms of NOs and NAs differ from their actual values for small n, the convergence is surprisingly fast when increasing n.
Although our calculations were done for the ground state only, the methodology can be applied mutatis mutandis to any other eigenstate (see [44] for their explicit forms) of this two-anyon system, resulting in the same NA asymptotic with suitably altered constant D(α).Our asymptotic results may also extend to eigenstates of anyon gases with higher numbers of particles, however proving this would require using a different set of mathematical tools such as the ones applied in the work [39].
The same method can be employed to derive the NOand NA-asymptotic for the non-interacting two-anyon system at hand in the fermion magnetic gauge.The fermionic counterpart of the wavefunction ( 8) is [42] , where θ 12 is the angle between z 1 and z 2 .Note that the coalescence cusp in Ψ F is of the same order as the coalescence cusp of Ψ B .By repeating the steps form Section III with Ψ B replaced by Ψ F we arrive at the identical conclusion concerning the asymptotics of NOs and NAs in the fermion magnetic gauge.
It is interesting to note that according to the power law (37), the anyonic NONs in 2D decay slower than for Coulombic multi-electron systems in 3D.This confirms the intuition that 2D anyon systems are characterised by strong correlations and, in light of work [40], are likely to be more challenging to tackle by the standard quantum chemistry toolset.
One might be tempted to interpret our presented results in terms of the Pauli exclusion principle for anyons.However, note that the magnetic gauge transformation (5) from the bosonic wavefunction Ψ B to the anyonyonic wavefunction Ψ α is nonlocal, thus the NONs of Ψ B are different from the NONs of Ψ α .Moreover, only the NONs of Ψ α are the ones that interpolate between fermionic case (Slater determinant for α = 1) and the bosonic case (fully condensed state for α = 0).Thus, for anyonic systems it is more appropriate to refer to gauge-invariant methods of measuring the Pauli exclusion principle such as the expectation value of the kinetic energy operator [56][57][58].

ACKNOWLEDGMENTS
The authors would like to thank Jonathan Robbins for helpful discussions.The numerical computations presented here were conducted using the University of Bristol HPC system (BlueCrystal 4 ).The research described in this publication has been funded by the National Science Center (Poland) under grant 2022/47/B/ST4/00002.The support from Max-Planck-Institut für Physik komplexer Systeme (MPI PKS), Dresden is also acknowledged by one of the authors (J.C.).
Appendix A: The leading-order expansion of Equation (25) The paper [48] provides the following theorem: for any f being an analytic function in the region Ω = {z ∈ C : a ≤ ℜ(z) ≤ b, ℑ(z) ≥ 0} we have

5 FIG. 2 .
FIG. 2. The NAs of the l = 0 sector for α ∈ 1 5 , 2 5 , 3 5 , 4 5 , 1 and M = 400.The number n represents the index of the values σn.The solid lines show the outcome of the fifth order polynomial regression for each value of α.

FIG. 4 . 1 / 5 .
FIG.4.Top row: the l = 0 natural orbitals of index 1 (a), and 6 (b) for various values of the fractional statistics parameter α.(c) shows a comparison of the numerical and asymptotic forms of the natural orbital n = 30, for which l = 0 and α = 1/5.The natural orbitals are observed to cross the axis times, where n is the principal quantum number of the orbital.Bottom row: a comparison of the numerically computed natural orbitals to their corresponding asymptotic formulae for large n, all with l = 0 and α = 1/5.For small r the deviations negligible.Frames (d), (e) and (f) show the oscillatory component of the natural orbitals, with n = 10 (d), n = 40 (e) and n = 80 (f).For small n, the deviations from the asymptotic formula diverge as r increases from 0. This divergence is less pronounced for larger n.