Exact closed-form analytic wave functions in two dimensions: Contact-interacting fermionic spinful ultracold atoms in a rapidly rotating trap

Exact two-dimensional analytic wave functions for an arbitrary number $N$ of contact-interacting lowest-Landau-level (LLL) spinful fermions are derived with the use of combined numerical and symbolic computational approaches via analysis of exact Hamiltonian numerical diagonalization data. Closed-form analytic expressions are presented for two families of zero-interaction-energy states at given total angular momentum and total spin $0 \leq S \leq N/2$ in the neighborhood of the $\nu=1$ filling, covering the range from the maximum density droplet to the first quasihole. Our theoretical predictions for higher-order spatial and momentum correlations reveal intrinsic polygonal, multi-ring crystalline-type structures, which can be tested with ultracold-atom experiments in rapidly rotating traps, simulating quantum Hall physics (including quantum LLL skyrmions).


II. METHODOLOGY
Extensions of Girardeau's mapping between impenetrable bosons and non-interacting spinless fermions [1], and similar mappings [8,10] applied to spinful and spinparallel fermions, led to the formulation of a hard-core boundary condition for strongly-repelling 1D fermions [11,12]. This entails vanishing of the many-body wave functions when two fermions with antiparallel spins are at the same position (in addition to the vanishing for parallel spins due to the Pauli exclusion principle). Concomitant of this condition is the appearance of 0IE states.
In CI calculations and for a given number N of spinful LLL fermions, the 0IE states emerge in each spin sector (S, S z ); see Fig. 1 for the case of N = 9 LLL fermions interacting with a i<j δ(z i − z j ) two-body potential, where z i = x i + iy i (with i = 1, 2, . . . , N ). Interest in such 0IE states arises from: (i) They can be prepared  (1). Spectra in a given spin sector S = 1/2, 3/2, 5/2, 7/2 and S = 9/2 are expilcitly denoted. The spectra were calculated for Sz = 1/2; however, they are independent of the precise value of Sz. The 0IE states of family A are colored in red; those of family B are colored in white. The first quasi-hole state is also explicitly denoted. L is the total angular momentum.
experimentally [39]. (For bosons, the experimental expectations for 0IE states include fittingly the ν = 1/2 bosonic Laughlin state [34,35,61].) (ii) They represent microscopic states that describe quantum LLL skyrmions [45,50]. (iii) The Laughlin states are 0IE eigenstates associated with short-range pseudopotential-type Hamiltonians [44,62]. (iv) For fully polarized fermions, 0IE states have been associated with the gapless edge excitations of the Laughlin droplet [63] in extended semiconductor samples. The many-body Hamiltonian describing ultracold neutral atoms in a rapidly rotating trap [29, 33-35, 39, 58] is given by where ω and Ω are, respectively, the parabolic trapping and rotational frequencies of the trap, and L denotes the total angular momentum, L = N i=1 l i , normal to the rotating-trap plane; the energies are in units of ω and the lengths in units of the oscillator length Λ = /(M ω), with M being the fermion mass. The first and second terms express the LLL kinetic energy, H K , and the third term represents the contact interaction, H int .
Our methodology integrating both numerical (e.g., fortran) and symbolic (algebraic, e.g., MATHEMATICA [64]) languages consists of two steps: (1) numerical diagonalization of the Hamiltonian matrix problem employing the ARPACK solver [65,66] of large-scale sparse eigen-value problems, followed by step (2) where the numerically exact CI wave functions are analyzed and processed using symbolic scripts targeting extraction of the corresponding exact analytical wave functions.
The basis Slater determinants that span the Hilbert space are where r, s = 1, . . . , N , the LLL single-particle orbitals are and σ signifies an up (α) or a down (β) spin. The master index I counts the number of ordered arrangements (lists) {j 1 , j 2 , . . . , j N } under the restriction that 1 ≤ j 1 < j 2 < . . . < j N ≤ K; K ∈ N is chosen large enough to provide numerical convergence. Below, explicit mention of the Gaussian factor is omitted.
Step (2) starts with the rewriting of the CI wave function Φ CI in Eq. (2) as where the replacement of the subscript "CI" by "alg" corresponds to the fact that, using the symbolic language code, one obtains an equivalent multivariate homogeneous polynomial Φ alg with algebraic coefficients c alg ; see the transcription of coefficients for N = 4 and N = 9 in Tables STI and STII in the SM [59].
Invariably, the symbolic code is able to simplify the derived multivariate polynomial in Eq. (5) to the compact form of a product of a Vandermonde determinant (VDdet), N i<j (z i − z j ), involving the space coordinates only, with a symmetric polynomial (under two-particle exchange) with mixed space and spin coordinates [see Eq. (6) below]. The factoring out of the VDdet reflects the fact that Φ alg represents a 0IE LLL state.
Using symbolic scripts, we verify further that the fullyalgebraic Φ alg [Eq. (6)] is indeed an eigenstate of the total spin, obeying the Fock condition [74]. The final closed form expressions [see Eq. (8) below] are derived for N ≤ 9, but they are valid for any N , thus circumventing the CI numerical diagonalization of large matrices, which is not feasible for N ≥∼ 10.
For the CI diagonalization, a small perturbing term V P (e.g., a small trap deformation [38], or a small hard-wall boundary [75]) needs to be added to the LLL Hamiltonian in Eq. (1). This has a negligible influence on the numerical eigenvalues, but it is instrumental in lifting the degeneracies among the 0IE states, and thus produce CI states whose total spin S is a good quantum number.

III. TARGETED TOTAL SPINS AND ANGULAR MOMENTA
For each size N , we provide analytic expressions for the maximum-spin (S = S z ) 0IE ground states with angular momenta L = L 0 + ∆L [with L 0 = N (N − 1)/2] from ∆L = 0 (maximum density droplet) to ∆L = N (first quasihole, 1QH); they form two families A and B (see Fig. 1 for an illustration).
Using k to denote the number of spin-up fermions and p that of spin-down fermions, and focusing on the case with k ≥ p (or equivalently p ≤ N/2), the states in both families are associated with the same set of total spins specified as S = S z = (k − p)/2 = N/2 − p. Furthermore, given a pair (k, p): (A) The states in family A have ∆L = p, with ∆L varying from 0 to N/2 for even N , and from 0 to (N −1)/2 for odd N .
(B) The states in family B have ∆L = k, with ∆L varying from N/2 to N for even N , and from (N + 1)/2 to N for odd N .
The states in family A are unique ground states, whereas those in family B are part of degenerate manifolds. (This degeneracy is lifted as described above.)

Mathematical preliminary
The quantity k-subset(list) is a subset containing exactly k elements out of the set of n elements (named list). The number of k-subsets on n elements is given by

General form of the 0IE LLL wave functions
The compact algebraic expression has the general form where χ(i) denotes an up spin, α, or down spin, β, and i = 1, . . . , N .
where l j = (j − 1) and i, j = 1, 2, . . . , N . The product of Jastrow factors above reflects the fact that the wave function in Eq. (6) is a 0IE eigenstate of the contactinteraction term, H int , in Eq. (1). Due to the fermionic symmetry of the Φ alg , Φ sym has to be symmetric under the exchange of any pair of indices i and j. Furthermore, Φ sym can be written as where P o m (defined below) are homogeneous multivariate polynomials of order o = p (family A) or o = k (family B), and  For each S = S z = (k − p)/2, except when k = p which has a single state, there exists a pair of targeted LLL states, with one state of the pair belonging to family A and the other to family B (see Fig. 1 for an example).
Family A: First, the following square matrices of rank p (the number of spin-down fermions) need to be considered: the matrices defined in Eq. (10). (Recall that k is the total number of spin-up fermions, and that {i 1 , . . . , i k } is also referred to as a second-level list.) The expression for the polynomial is given by where the symbol "Perm" denotes a Permanent.
The analytic expressions of the states with S z < S, in a given spin multiplicity 2S +1, are obtained by repeated application of the spin lowering operator.
Note that the first quasi-hole state (1QH) [41,75] coincides with the analytic expression associated with family B above for L = L 0 + N .

VI. CONCLUSION
A novel approach for deriving exact closed-form analytic expressions for the wave functions (beyond the Jastrow-factors paradigm) of an assembly of 2D contactinteracting spinful LLL fermions (for any N ) was introduced and validated. Such expressions require as input only the three parameters N (number of particles), L (total angular momentum), and S (total spin). Examples were presented for two families of zero-interaction-energy states, from the maximum density droplet to the first quasihole in the neighborhood of ν = 1. Ensuing theoretical predictions for higher-order momentum correlations for N = 19, 25, and 27, revealing intrinsic polygonal, multi-ring crystalline configurations, could be tested with ultracold-atom experiments in rotating traps simulating spinful quantum Hall physics, including LLL skyrmions. The present approach can be extended to the neighborhood of any ν = 1/m that starts with a Laughlin wave function.

VII. ACKNOWLEDGMENTS
This work has been supported by a grant from the Air Force Office of Scientific Research (AFSOR) under Award No. FA9550-21-1-0198. Calculations were carried out at the GATECH Center for Computational Materials Science.
Appendix: Comparison with the symmetric polynomials for quantum skyrmions in Ref. [45] We compare here with the symmetric polynomials for the seed skyrmions specified in Eq. (6) of Ref. [45] or Eq. (2) in Ref. [49].
Omitting the trivial Gaussian functions, these polynomials are given by the single formula where Z m are the spin primitives defined in Eq. As is the case in Ref. [45], one can take the index m as running over the p-subsets associated with the spin-down fermions, because there is a one-to-one correspondence to the k-subsets of the spin-up fermions. Note that Ref. [45] (Ref. [49]) uses the capital letter K (R) in place of our p. We consider the case N = 5, k = 4, p = 1, S = S z = 3/2, and ∆L = 4, belonging to family B in our exposition.
According to Eq. (A.1), the corresponding MFB symmetric polynomial is The corresponding exact symmetric polynomial derived in this paper is given by Eqs. (8) and (13)  Omitting the trivial Gaussian functions, these polynomials are given by the single formula where Z m are the spin primitives defined in Eq. (9) of the main text, i.e., The index m runs also over the k-subsets We present here a comparison for the case N = 5, k = 4, p = 1, S = S z = 3/2, and ∆L = 4, a state belonging to family B in our exposition.
According to Eq. (S5), the corresponding MacDonald-Fertig-Brey (MFB) symmetric polynomial is The corresponding polynomial derived in this paper is given by Eqs. (8) and (13) in the main text. Expanding the permanent, one obtains for the polynomial with k = 4 and m = 1 (in front of the Z 1 spin primitive): For the P 4 2 polynomial in front of Z 2 , one obtains similarly: In general, one has with P 4 m [z] = c 1 z 1 z 2 z 3 z 4 + c 2 z 1 z 2 z 3 z 5 + c 3 z 1 z 2 z 4 z 5 + c 4 z 1 z 3 z 4 z 5 + c 5 z 2 z 3 z 4 z 5 , and c i = 4 when i = m and c i = −1 otherwise. We note that the expressions associated with the Z i , i = 1, . . . , 5 in the MFB polynomial consist only of a single term with a numerical factor +1 in front. This contrasts with our expressions in Eqs. (S8), (S9), and (S11) which have five terms each with factors of +4 and -1 in front of them.
For p (spin-down fermions) > k (spin-up fermions), the MFB expression in Eq. (S5) is associated with a negative total-spin projection S z = (k − p)/2 < 0. In this case, the indices for the corresponding wave function in this paper are found by reversing all N spins, i.e., by considering the case with p → k, k → p, and S z = |(k − p)/2|.
Using our algebraic scripts, we readily verified that the wave functions derived in this work are indeed eigenfunctions of the total-spin square operator [with eigenvalue 15/4 and S = 3/2 for the case in this section], whereas the MFB ones are not (as indeed has been discussed by M. Abolfath et al., Phys. Rev. B 56, 6795 (1997) (https://doi.org/ 10.1103/PhysRevB.56.6795).
In particular, applying the spin-square,Ŝ 2 , operator, one gets On the contrary, for the MFB wave function, one gets Ŝ 2 − 15 4 Φ sk,MFB p=1 = (z 1 z 2 z 3 z 4 + z 1 z 2 z 3 z 5 + z 1 z 2 z 4 z 5 + z 1 z 3 z 4 z 5 + z 2 z 3 z 4 z 5 ) TABLE STI. The 15 dominant numerical CI coefficients, cCI(J), in the CI expansion of the relative LLL ground state, and the corresponding extracted algebraic ones, c alg (J), for N = 4 fermions with total angular momentum L = 8 (first 0IE state) in the (S = Sz = 0) spin sector of singlet states. The spinful-fermions Slater determinants are specified through the set of single-particle angular momenta and spins, (l1 ↑, l2 ↑, l3 ↓, l4 ↓). The converged CI expansion included a much larger set of basis Slater determinants, but naturally the dominant ones are only relevant, namely those with coefficients |cCI(I)| > 0.001.