Reduced density matrix approach to ultracold few-fermion systems in one dimension

The variational determination of the two-fermion reduced density matrix is described for harmonically trapped, ultracold few-fermion systems in one dimension with equal spin populations. This is accomplished by formulating the problem as a semi-definite program, with the two-fermion reduced density matrix being subject to well-known $N$-representability conditions. The ground-state energies, as well as the density, pair-correlation function, and lower-order eigenvalues of the two-fermion reduced density matrix of various fermionic systems are found by utilising an augmented Lagrangian method for semi-definite programming. The ground-state energies are found to match well to those determined by full-configuration interaction and coupled-cluster calculations and the density, pair-correlation function, and eigenvalue results demonstrate that the salient features of these systems are well-described by this method. These results collectively demonstrate the utility of the reduced density matrix method firstly in describing strong correlation arising from short-range interactions, suggesting that the well-known N-representability conditions are sufficient to model ultracold fermionic systems, and secondly in illustrating the prospect of treating larger systems currently out of the reach of established methods.


I. INTRODUCTION
The study of few-fermion systems in quantum mechanics began in earnest immediately after the discovery of the Schrödinger equation in 1925 [1]. The first systems to be analysed were simple atoms and molecules. The treatments were extended to larger systems as the computational power increased and the theoretical methods became more well-founded. In more recent years, experimental data have permitted the study of more exotic few-fermion systems such as quantum dots [2][3][4][5] and trapped, highly-correlated ultracold fermionic gases [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. The ongoing theoretical efforts to accurately and efficiently capture the ground-state properties and dynamics of these systems coupled with ever-growing experimental capacities has ensured that the field of few-fermion systems retains immense interest.
A main goal of many-body quantum mechanics is the accurate description of correlation in a system, which constitutes the intricate interactions that are not captured by a mean-field approach. In typical atomic and molecular systems, long-range Coulombic interactions are responsible for correlation in the system. However, in ultracold systems, short-range, point-like interactions dominate and it is these characteristic interactions that give rise to the exotic effects exhibited by these systems. An apt many-body scheme would naturally capture the effects of correlation in a system regardless of the nature of the interaction. However, there is no a priori reason that a scheme which accurately captures correlation from long-range effects will also capture correlation from short-range effects to the same degree of accuracy.
Fortuitously, many schemes that arose out of a need to capture correlation arising from Coulombic interactions are equally applicable to short-range interactions. Many * knightm1@student.unimelb. edu.au of these schemes are known as post-Hartree-Fock (HF) methods as they utilise the non-correlated HF groundstate as an initial step to construct a wavefunction solution to the Schrödinger equation. For example, the coupled-cluster (CC) method has seen wide success in atomic and molecular calculations [26][27][28][29][30][31][32] and, more recently, has been shown to accurately capture the groundstate of ultracold bosonic [33,34] and fermionc systems [14,35], thus demonstrating the practical generality of the method.
Included in the myriad techniques that constitute post-HF methods is the technique of full-configuration interaction (FCI) which variationally finds an exact solution of the Schrödinger equation in the Hilbert space spanned by the finite basis set used. This method is equivalent to a direct diagonalisation of the Hamiltonian and serves as a benchmark against which the accuracy of other methods can be compared. Naturally, it has been applied extensively in quantum chemistry in particular and many-body physics more generally. However, it must be noted that such an approach is computationally expensive, and cannot be extended to arbitrarily large systems.
Due to the immense number of parameters in the Nfermion wavefunction, for moderate to large N in FCI calculations, where N is the number of fermions in the system, it is enticing to consider whether the groundstate properties of these systems can be found using reduced quantities. These quantities would ideally depend on fewer parameters than the N -fermion wavefunction. Density functional theory (DFT) is one such approach, which utilises the single-particle electron density as the fundamental variable rather than the N -fermion wavefunction. It has seen wide success in capturing correlation in atomic and molecular systems, which is the correlation that arises from the Coulombic interactions, as well as larger, more complicated systems such as those emerging in materials science [36][37][38][39]. DFT also benefits from excellent scalability, permitting its application to large-scale systems as, generally speaking, the approach does not scale with particle number, N . The reduction in computational complexity attained by utilising the electron density, or, more generally, the fermion density, comes at a cost, however, as the exact energy functional for the system is replaced by an unknown energy functional of the density. This replacement brings about numerous issues concerning the accuracy and range of applicability of DFT [40].
An alternative approach to DFT which still uses reduced quantities utilises the two-fermion reduced density matrix (2-RDM), which may be derived from the N -particle density matrix by integrating out the degrees of freedom of all but two of the fermions. The utility of the approach arises from the property, first demonstrated by Husimi [41], Löwdin [42], and Mayer [43], that any quantum system of N particles with at-most pairwise interactions can be described completely by the 2-RDM [41]. In such a case, the variational object is the 2-RDM and the energy functional to minimise is the Hamiltonian itself. In other words, physical characteristics of the system such as the ground-state energy, particle density, and the pair-correlation function can be expressed as linear functionals of the 2-RDM. Earlier, Dirac [44] had shown that the HF ground-state for an N -fermion system could be expressed explicitly in terms of the one-fermion reduced density matrix (1-RDM). Soon thereafter, the idea of reducing an N -fermion problem to an effective two-particle problem was promoted by A. J. Coleman [45]. Early attempts by Mayer [43], Tredgold [46], and Coleman [45,47] to utilise this methodology yielded poor results, with the ground-state energy being far too low. Indeed, it was soon realised that a Rayleigh-Ritz minimisation of the ground-state energy subject to variation in the elements of the 2-RDM failed to impose sufficient constraints.
These additional constraints were termed Nrepresentability conditions [47], and were required to ensure that a trial 2-RDM that is determined variationally corresponds to a legitimate N -fermion wavefunction. The search for a complete set of such conditions was known as the N -representability problem and was pursed actively in electronic structure theory and many-body quantum mechanics for more than five decades. While the complete solution to the N -representability problem eluded the theoretical community, Coleman [47] and Garrod and Percus [48] obtained certain constraints known as the D, Q, and G conditions. These manifested as semidefinite constraints on matrices representing the probability distribution of two fermions, one fermion and one hole, and two holes, respectively. While the D, Q, and Q conditions are necessary to constrain the 2-RDM to be N -representable they are only sufficient in the N = 2 fermion case and additional conditions are needed for larger systems. Further conditions were soon found by Erdahl [49], which are known as the T1 and T2 conditions, which constitute semidefinite constraints on matrices corresponding to the various distributions of three fermions and holes.
Motivated by the development of sophisticated and efficient methods in semidefinite programming (SDP) [50,51], significant progress was made during the early 2000s to apply these N -representability conditions while variationally determining 2-RDMs. The development of powerful interior-point methods in the late 1990's spawned the rise of interest in the 2-RDM method as a whole and, later on, significantly more powerful boundary-point methods were developed, which greatly increased the efficiency of the implementation [52][53][54]. Examples of systems that have been analysed using this formalism include simple atomic and molecular systems [52,53,[55][56][57][58][59][60] and quantum dots [61] in which correlation primarily arises from the long-range Coulombic interactions in the system and also the Hubbard model [60,[62][63][64], in which short-range interactions are present. Moreover, the RDM method has been recently applied to a large range of transition metal systems [65,66] as well as strongly correlated molecules arising in both organic and inorganic chemistry [67][68][69].
In this work, we apply the RDM method to the field of trapped, ultracold few-fermion systems exhibiting a contact interaction in both the weakly-and stronglyinteracting regimes. For simplicity, we consider onedimensional systems of net zero spin yet we permit the interactions in these systems to be repulsive and attractive. Such a system serves as an excellent medium in which to demonstrate the effectiveness of the RDM method in capturing short-range interactions. In particular, the analysis of short-range, point-like interactions with the RDM methodology permits us to assess the accuracy with which correlation in these systems is captured by the commonly utilised N -representability conditions which, until now, have predominately been utilised in the study of long-range, Coulombic interactions. Moreover, the strength of the contact interaction in these ultracold systems can be made arbitrarily large, thus we can consequently analyse the effectiveness of the RDM method in capturing substantial correlation in these systems in both the attractive and repulsive regimes.
Ultracold quantum gases are of immense interest as they serve as a testbed for the study of many-body quantum mechanics more generally, as they are highly controllable [17,70]. Also, the intricacies of the interactions in simple ultracold gases with a small number of particles serve primarily to inform the ground-state properties and dynamics of larger, and more complicated condensed matter systems [71]. Thus, the application of quantum chemistry techniques to trapped ultracold fewfermion systems has allowed the treatment of larger systems where an exact diagonalisation of the Hamiltonian is completely unfeasible [14,35]. In this study, the utility and accuracy of the RDM method for highly-correlated systems is demonstrated in the ground-state calculations for simple N = 2, 4, 8, 10, and 12 fermion systems.
The organisation of this paper is as follows. In Sec. II we present the general theory of the RDM approach, and demonstrate how the expectation value of a many-body Hamiltonian can be expressed directly in terms of the 2-RDM. In Sec. III we briefly outline the N -representability problem, and discuss the constructive solution discovered by Mazziotti [54] while explicitly giving the matrix representations of the well-known conditions. In Sec. IV we describe the implementation of the RDM as an SDP. In Sec. V we briefly discuss ultracold few-fermion systems and in Sec. VI we demonstrate the results of the RDM calculations. We conclude with a discussion on the limitations of the RDM methodology, as well as discussing the key benefits of using the RDM in the ultracold fewfermion field of research.

II. REDUCED DENSITY MATRIX THEORY
Consider an N -fermion quantum system. It is an axiomatic assumption that the Hamiltonian for this system can be expressed as the sum of one-body and two-body operators as, where (1)ĥ (i) is a one-body operator that acts on a single fermion, i, and (2)V (i, j) is a two-body operator that acts on the pair of fermions, i and j.
Let |Ψ be the (pure) state vector describing the ground-state of this system. Then, the ground-state energy of the system, E 0 , is determined from the timeindependent Schrödinger equation [1],Ĥ |Ψ = E 0 |Ψ . The N -fermion density operator, (N )ρ , for this system, as introduced by von Neumann [72], is given by the outer product of the state vector |Ψ with itself, (N )ρ = |Ψ Ψ|. Projecting this operator into coordinate space yields the N -fermion density matrix, where x i denotes the spin and spatial coordinates of the i th fermion [73].
The diagonal element of the N -fermion density matrix, when x 1 = x 1 , . . . , x N = x N , has the usual statistical interpretation: (N ) D(x 1 , . . . , x N ; x 1 , . . . , x N ) is a probability distribution function describing the locations in space and spins of the N fermions.
The diagonal element of a p-fermion reduced density matrix (p-RDM) is the probability distribution function for p < N fermions occupying any p of the N positions and spins given by x 1 , . . . , x N , and is found by integrating out the spin and spatial coordinates of the (p + 1)th to the N th fermions, This object was first introduced by Husimi [41] although Dirac [44] showed in 1930 that a Hartree-Fock groundstate for a fermionic system can be expressed in terms of 1-RDMs. The normalisation factor p! N p was introduced by McWeeny [74]. Other normalisation factors exist in the literature today, with the factor N p due to Löwdin [42] and the factor of unity due to ter Haar [75] and Coleman [47].
A key advantage in employing p-RDMs lies in a fundamental result that the expectation value of a p-body operator can be expressed as a linear functional of the p-RDM [41,42,47,74,75]. Since the Hamiltonian in Eq. (1) is the sum of one-and two-body operators, the ground-state energy, E 0 = Ψ|Ĥ|Ψ , can be written as a linear functional of the 1-and 2-RDMs, which are given by To apply the RDM apparatus to physical problems, a finite basis set must be introduced. Typically, this is a single-particle basis of spin-orbitals, which we denote by {ϕ i (x)}. Let us denote the size of the spin-orbital basis set by K. Due to the two-fold spin multiplicity of fermions, a rank K spin-orbital basis means we have K/2 spatial-orbital functions in the given basis. In such a basis, the 1-and 2-RDMs are represented by In Eq. (6) and (7) the tensors (1) D i j and (2) D ij kl can be expressed conveniently in second-quantised notation, providing the most common expressions for the 1-and 2-RDMs as found in the literature, whereâ † andâ are the annihilation and creation operators for fermions. We also note that we can express the general Hamiltonian in Eq. (1) in a second-quantised form aŝ where the indices i, j, k and l correspond to elements of the spin-orbital basis set {ϕ i (x)}, the sums run over all possible values of the indices in the basis, i, j, k, l ∈ {1, 2, . . . , K}, and (1)ĥi j and (2)V ij kl are the one-and twofermion integrals, (1) By taking the expectation value of the Hamiltonian in Eq. (10) with respect to the ground-state |Ψ and incorporating the expressions of the 1-and 2-RDMs in secondquantised notation (Eqs. (8) and (9)) we see that the ground-state energy is a linear functional of the 1-and 2-RDMs;

III. N -REPRESENTABILITY
Initial efforts to calculate the ground-state properties of simple systems with the RDM method yielded drastically incorrect results [45][46][47][48]76]. This is due to the lack of N -representability conditions applied to the 2-RDM. The search for these N -representability conditions, which ensure that the 1-and 2-RDMs correspond to a legitimate N -particle wavefunction, are what Coleman [47] called the N -representability problem.
We denote by (N ) P the set of all positive semidefinite Hermitian operators of unit trace on the Hilbert space for our N -fermion system. This set is equivalent to the set of all possible N -fermion density operators (N ) D [47,77]. We also denote by (p) P N the set of all N -representable p-RDMs. The prescription given in Eq. (3) permits one to calculate a p-RDM from the N -particle density matrix but does not provide a mechanism to ensure that an alleged p-RDM corresponds to a legitimate N -particle density matrix. That is, a variationally determined p-RDM does not necessarily belong to the set (p) P N but to the more general set of p-RDMs which may not be Nrepresentable, the set of which we denote (p) P. Hence, (p) P N ⊂ (p) P, i.e., (p) P N is a proper subset of (p) P. Thus, the N -representability problem can be stated as the search for a complete characterisation of the set (p) P N as a subset of (p) P. Now, the set (N ) P is convex as are its subsets (p) P and (p) P N , with (p) P N being a convex subset of (p) P. A convex set is a set S such that if x, y ∈ S then αx + βy ∈ S, for all α, β ∈ R ≥0 such that α + β = 1. That is, convex combinations of density matrices are themselves density matrices. An extreme element of a convex set is an element that cannot be expressed as a convex combination of other elements in the set [78,79]. These are the density matrices representing pure states as opposed to mixed states, which are themselves represented by convex combinations of pure states. By the Krein-Milman theorem [80], a convex set is completely specified by its extreme elements. From this theorem, a complete characterisation of the set (p) P N can be garnered from its extreme elements, i.e., by considering the pure state density matrices of the system.
The complete characterisation was first discovered by Mazziotti [81,82]. Define the convex set of p-fermion operators that are positive semidefinite in their trace against an N -representable p-RDM by (p) P * N , Kummer [83] showed that this set explicitly exists, and this set provides a complete characterisation of its dual or polar set, i.e., the set of N -representable p-RDMs, (p) P N . As such, a knowledge of the extreme elements in (p) P * N allows us to completely characterise the set of p-RDMs.
Mazziotti noted the important relation that (2) where K is the rank of the spin orbital basis. As such, the extreme elements of (2) P * N are specified by convex combinations of the extreme elements of (K) P * N , which are of the form i is a polynomial in creation and annihilation operators of order ≤ K. Hence, the operators which constrain the 2-RDM to be N -representable, for example, are represented by (2)Ô = i ω iĈiĈ † i where the weights ω i are chosen such that all three-body or higher operators cancel on their trace against the 2-RDM.
Such a solution yields a collection of (p, q)-positivity conditions [54], where the p indicates the highest RDM required to evaluate the condition (here that remains the 1-, and 2-RDM) and the q indicates the highest qbody operator cancelled by the combinations expressed above. For example, the (1, 1)-positivity conditions are the conditions on the 1-RDM which come from consider- Keeping the trace of these operators with the ground-state wavefunction greater than or equal to zero yields the Nrepresentability conditions on the 1-RDM, first found by Coleman [47]; where 0 means 'positive semiefinite', δ is the Kronecker delta, and Eq. (15) follows from the anticommutation relation for fermion creation and annihilation operators.
The constraints here are necessary and complete, that is, the 1-RDM is completely N -representable if equipped with these constraints, and the constraints are equivalent to the Pauli exclusion principle: the occupation numbers of a given state in a fermion system must be 0 or 1 [47]. The (2, 2)-positivity conditions arise from considering (2)Ô G then gives us the wellknown D, Q, and G conditions, first due to Coleman [47] and Garrod and Percus [48], The conditions given in Eq. (16), Eq. (17), and Eq. (18) are physically interpreted as constraining the probability distributions of two fermions, two holes, and one fermion and one hole to be non-negative, respectively. They do not form a complete set, however, of Nrepresentability conditions on the 2-RDM. For example, the (2, 3)-positivity conditions are also necessary for a description of a system of N > 2 fermions. Combinations of constraint matrices derived from consideringĈ i operators of degree 3 are utilised, since these can be expressed as linear combinations of the elements of the 1and 2-RDMs. These are called the T1 and T2 conditions, and have the expression [49], where These four matrices corresponding to the distributions of three fermions, two fermions and a hole, one fermion and two holes, and three holes, respectively. Although these matrices are created through consideration of threebody operators, the combinations given by Eq. (19) and Eq. (20) constrain the 2-RDM through cancellation. The T1 and T2 conditions were found implicitly by Erdahl [49] and first implemented by Zhao [58] and Mazziotti [56,57]. All conditions expressed here can be recast as linear constraints on the 1-and 2-RDMs by implementing the fermion anti-commutation relations. For example, the (2) Q and (2) G matrices can be expressed as, Similar expressions can be found for the (3) T 1 and (3) T 2 matrices are given elsewhere (for example, see Zhao et al. [58]). With these expressions, a set of linear constraints on the components of the 1-and 2-RDMs are now available in a variational calculation. In this work, we include conditions up to and including the T1 and T2 conditions and therefore, except in the N = 2 fermion case, we hold the 2-RDM to be be approximately N -representable. We will see that the (2, 2)and (2, 3)-positivity conditions, given by constraining (2) D, (2) Q, and (2) G, and (3) T 1 and (3) T 2 to be positive semidefinite, yield accurate results for the ground-state energy at a range of short-range interaction strengths when compared to other methods.

III.1. Spin-Adaption
The basis set we introduced is such that each spinorbital is a product of a spatial wavefunction and a number representing the spin, spin-up or spin-down, respectively, where α and β are eigenfunctions of the spin operatorŜ z . That is, we have L spatial wavefunctions and hence K = 2L spin-orbitals, allowing for the twofold spin multiplicity. We also order our basis appropriately such that the first L spin-orbitals all are spin-up (or all spin-down) and the next L are all spin-down (or all spinup). Then, the spin-orbitals indicated by the index 1 and L + 1 have the same spatial wavefunction, but differ in their spin.
Through this construction of the basis, we render our 2-RDM and associated matrices block-diagonal in a process known as spin-adaption. Typically, the term spinadaption refers to constraining the basis functions to be simultaneous eigenfunctions ofŜ z andŜ 2 [84], but here we simply enforce the former condition and implement the latter as a constraint enforced during the variational optimisation through Eq. (37). This form of implementation has been used numerous times in applications to atomic and molecular systems [55, 56, 58-60, 85]. Spin and symmetry adaption from both theŜ z andŜ 2 operator of the (2) D, (2) Q, and (2) G matrices has been analysed [86].
The spin-adapted 1-RDM has the form, Here, the blocks are identical due to the ordering of the spin-orbital basis, and have size (K/2) × (K/2). The spin-adapted 2-RDM has the form, In the case of a closed-shell system, where each spatial orbital is doubly-occupied, the first and second blocks of the 1-and 2-RDMs are identical, and the sizes of the first two blocks of the 2-RDM are K(K/2 − 1)/4 × K(K/2 − 1)/4, whereas the third block has size K 2 /4. The (2) Q matrix has the same block structure as the 2-RDM since it is the two-hole reduced density matrix and holes obey the same statistics. The one-hole one-fermion RDM, (2) G, has a more complicated block structure as it does not exhibit as many symmetry properties as the other matrices. It has the spin-adapted form, The block structure of the (3) T 1 and (3) T 2 matrices have more complicated representation which we will not give here (see Zhao et al, [58]).

III.2. A Summary of the N -Representability Conditions
In addition to the N -representability conditions which manifest in the D, Q, G, T1, and T2 conditions, we have the conditions that hold for density matrices in general. We list these conditions here.
i) The 2-RDM is anti-symmetric in its upper and lower indices, The (2) Q matrix has exactly the same antisymmetry as the 2-RDM, whereas the (2) G matrix does not. Furthermore, the (3) T 1 matrix is antisymmetric with respect to interchange of its triples of indices, and the (3) T 2 matrix is anti-symmetric with interchange of the latter pairs in each triple.
ii) Hermiticity of the 1-and 2-RDMs, This Hermiticity is also extended to the (2) Q, (2) iii) We ensure that the constant number of fermions (spin-up and spin-down) in the physical system is reflected in the trace of the 1-and 2-RDMs. That is, We recall that this normalisation scheme is a convention, as the 1-and 2-RDMs can be normalised to unity.
iv) The number of spin-up (or spin-down) fermions is held constant, a constraint which must show up in the normalisation: where i α , j α , . . . indicate indices corresponding to spin-up orbitals, and N α is the number of spinup orbitals. The latter of these constraints can be derived from the fact that the wavefunction, |Ψ , of the system is an eigenstate of the number operator for the spin-up electrons,N α [55].
v) We ensure we have a linear constraint involving the net spin, S. This is due to the fact that |Ψ is held to be an eigenfunction of theŜ 2 operator. This condition is difficult to express in the notation above, so we use an alternate notation to ensure this constraint is explicit. Following Zhao et al. [58], we let n i , n j , . . . denote spatial orbital indices, with n i α and n i β denoting spin-up and spin-down spatial orbitals, respectively. Then, we ensure ni,nj Note that in all applications in this work, the net spin is always zero, i.e., S = 0.
vi) We ensure we have the correct spinsymmetries of the matrices where σ i is the spin of the ith spin-orbital. There are similar constraints for the remaining matrices, which are listed elsewhere [58].
vii) The N -representability constraints; the posi- Due to the symmetries (and anti-symmetries) of the Hamiltonian and 1-and 2-RDMs, as well as the other matrices that are constrained to be positive semidefinite, we can restrict the elements of the matrices we use in the calculation to increase efficiency. For example, due to the anti-symmetry of the 2-RDM, we can now only consider elements of (2) D ij kl for 1 ≤ i < j ≤ K and 1 ≤ k < l ≤ K. We can do the same for (2) Q ij kl . However, the (2) G matrix does not permit such a reduction as it does not have the same anti-symmetric structure. For (3) T 1 ijk lmn we can take 1 ≤ i < j < k ≤ K and 1 ≤ l < m < n ≤ K, and for (3) T 2 ijk lmn we take those elements with 1 ≤ j < k ≤ K and 1 ≤ n < m ≤ K, where K is the spin-orbital basis rank.
Similarly, in the Hamiltonian given in Eq. (10), with the matrix elements of the two-body operator, (2) V ij kl , we take the elements 1 ≤ i < j ≤ K and 1 ≤ k < l ≤ K. As such, the optimisation problem becomes, where (2)Ṽ is now the anti-symmetrised two-fermion integral,

IV. IMPLEMENTATION AS A SEMIDEFINITE PROGRAM
The type of optimisation problem given by Eq. (39) with the necessity to constrain the constituent matrices to be positive semidefinite is one that naturally facilitates the use of semidefinite programming. A semidefinite program (SDP) is a general optimisation problem that can be expressed as [50,51], where h, y ∈ R m are vectors (where m is the number of constraints), and A p , C ∈ S n are real symmetric matrices of size n (where n depends on the number of constraints and the basis rank. See Fukuda et al, [59]). Here, h T denotes the transpose of h, such that the product h T y represents the typical inner product between two Euclidean vectors.
We formulate the RDM variational problem as an SDP by letting the vector h contain the matrix elements of the Hamiltonian and the vector y contain the components of the 1-and 2-RDM which will be determined variationally. Then, by judicious choices of the matrices A p and C we can implement the conditions (1)-(7) given in the previous section. Explicitly, let us define a linear transformation: svec : S n → R n(n+1)/2 , where S n is the set of real, symmetric matrices of size n × n, by where U ∈ S n . The appearance of the √ 2 ensures that we need not consider two products of the form, for example, (1) D 2 1 , but simply double-count one of them as they are equal. This transformation allows us to consistently place the elements of a symmetric matrix into a vectorised form. Hence, in the SDP above, we let h = (svec( (1) h i j ) T , svec( (2)Ṽ ij kl ) T ) T and similarly for the y vector which contains elements of the 1-and 2-RDMs [58,59]. We then ensure that the matrix held to be positive semidefinite, which is constructed by the A p and C matrices, has the following diagonal blocks: Hence, if we consider the first block in which we are constraining (1) D to be positive semidefinite, we note that the first elements of the y matrix are those containing the elements of (1) D, i.e., y = ( (1) D 1 1 , √ 2 (1) D 1 2 , . . . ) then one can see that in the corresponding blocks of A p , we would have

  
and the first block of the C matrix would contain all zeros. As such, the matrices involved in such a calculation are extremely sparse. The size of sparsity of the SDPs for RDM implementation has been covered in detail elsewhere [59].
The formulation of the RDM apparatus as an SDP was implicit in the early works of Garrod et al. [87,88] on atomic beryllium and Rosina et al. [89]. The development of powerful interior points methods [90,91] for SDPs were developed largely in the late 1990s which led to a resurgence in the application of the RDM theory to atoms and molecules [55,56,58]. These calculations yielded accurate results, but were restricted to relatively small basis sizes. However, boundary-point methods [92] have facilitated an enormous increase in computational efficiency, allowing the analysis of more complicated systems [52][53][54]. In this work, we utilise the Semidefinite Programming Algorithm (SDPA) [93,94] which utilises interior point methods, as well as SDPNAL+ [95][96][97] which utilises an augmented Lagrangian method [95,98], rendering it extremely efficient for large-scale SDPs. To implement these programs, we utilised Spartan highperformance computing system at the University of Melbourne [99].

V. APPLICATION TO ULTRACOLD FEW-FERMION SYSTEMS
In this work, the RDM method is not spruiked as a method leading to extraordinary accuracy relative to available methods in the field, but a promising alternative to treat large systems that is competitive in its accuracy and efficiency to known methods whatever the basis used. Due to the method's dependence on the basis rank and not the particle number N , improvements in optimisation will allow small to moderately sized systems that are highly correlated to be analysed.
The utility of the RDM method has been demonstrated in the field of quantum chemistry but has seen few applications outside this field. Here, we apply the method in another field, that of ultracold few-fermion systems. We demonstrate that the RDM method with no more constraints than the T1 and T2, which are widely implemented in quantum chemistry, are sufficient to capture strongly correlated ultracold systems exhibiting contact interactions even as the system approaches the unitary regime, where the interaction strength becomes arbitrarily large.
We choose the most simple of such systems: a onedimensional, harmonically trapped system of N particles with even populations of spin-up and spin-down particles with a short-range s-wave interaction. These systems serve as a convenient medium to assess the accuracy of the RDM method, and the utilisation of the T1 and T2 conditions, in capturing strong correlation in general ultracold systems, as the interaction strength can be made arbitrarily large, and the intrinsic nature of the point-like contact interaction exhibited in this system is not peculiar to one-dimensional systems, permitting us to expect comparable accuracy in other dimensions. Such systems have been studied at length both theoretically [6-19, 35, 100] and experimentally [19][20][21][22][23][24][25]. Naturally, quasi-one-dimensional systems are constructed by the means of relatively strong harmonic confinement in two orthogonal directions and a relatively weak confinement in a third. If the confinement in the planar direction(s) is strong enough, we can model the system with a one-dimensional Hamiltonian.
By a comparison to Eq. (1), we make the identifica- Here, we are working in units where = m = ω = 1 where m is the mass of the fermion and ω is the onedimensional trapping frequency. Also, g 1D parametrises the interaction strength and has the units ωa ho where a ho is the oscillator length of the harmonic trapping potential. The parameter g 1D is related to an effective one-dimensional scattering length a 1D via g 1D = −2 2 /(ma 1D ). This one-dimensional scattering length can be derived from the s-wave scattering length, a s , and the characteristic length of the planar trapping potential a ⊥ through the equation where C = 1.4603 [101]. Thus, while the parameter g 1D serves as a measure of the strength of the contact interactions in our system, it is also dependent on the s-wave scattering length, a s , and the trap geometry. The energy of this system is consequently measured in units of ω.
A convenient basis to study this system is the set of one-dimensional quantum harmonic oscillator eigenfunctions which are defined by, where H n (x) is a Hermite polynomial. The error in many-body calculations using this basis scales approximately as 1/ (K/2) [14] where we recall K is the spinorbital basis rank. Other basis sets exhibit superior scaling, such as the plane-wave basis which scales as 1/(K/2) [102]. However, we utilise the single-particle wavefunctions in Eq. (45) due to the simplicity of the implementation (in the harmonic oscillator basis, the matrix elements of the Hamiltonian can be calculated quickly, as the one-electron integrals are given by the well-known harmonic oscillator energy i| (1)ĥ |j = δ ij (i + 1/2), and the two-electron integrals can be calculated exactly using Gauss-Hermite quadrature) and to allow this paper to be easily compared to works with similar result utilising the same basis (see, for example, [13][14][15]35]). With this choice of basis, the results presented here do not exhibit any improvement in the accuracy of the ground-state energies of N fermion systems relative to results already presented in the literature, yet the results do show that the RDM method with appropriate Nrepresentability conditions accurately captures the correlation in these systems when compared to other wellknown, and accurate, methods.
To assess the accuracy of the ground-state energy found from the RDM method, we compare it to the exact ground-state energy of the system, in the given basis, if feasible. This is found from the FCI method, which is equivalent to a direct diagonalisation of the Hamiltonian. We note that configuration-interaction methodologies peculiar to ultracold systems with contact interactions have been developed [103,104], yet we will compare to a more general FCI method as the latter is more often utilised and is more general in scope.
FCI calculations are infeasible for large systems, since an upper bound to the scaling of the computational complexity of the FCI problem is (K!/(N/2)! (K − N/2)!) 2 [105]. For such systems where the FCI result is not available, we compare to the energy as found from coupledcluster (CC) method [26][27][28][29][30][31][32]. This method utilises the HF ground-state and accounts for correlation through the use of a cluster operator, which acts on the reference wavefunction and produces linear combinations of excitations [32]. The inclusion of single and double excitations is called the CCSD method, and if we further include approximate contributions from triple excitations which are found using many-body perturbation theory, we have the CCSD(T) method, which can be regarded as the gold standard in terms of efficiency and accuracy for small to medium size systems. Note that the CC method has been used for the one-dimensional ultracold few-fermion system extensively in the work by Grining et al. [14,35]. In this work, the FCI, CCSD, and CCSD(T) results were found by using the PySCF program, a Python-based quantum chemistry package [106] where we customised the Hamiltonian for our purposes.
Because of the highly-degenerate nature of the larger systems considered in this work the ground-state energy itself does not serve as a conclusive metric of the accuracy of the RDM energy at larger interaction strengths nor for larger systems [107][108][109]. As such, a comparison to known results in the field in such a case is pertinent. In the presence of infinitely strong repulsive interactions, i.e., when g 1D → +∞, the N -fermion wavefunction can be derived exactly [7,19,[107][108][109]. This is because the infinitely strong repulsive interaction between two fermions can be cast as a manifestation of the Pauli exclusion principle, except applying to fermions without the same quantum numbers. Then, any wavefunction of the form Ψ(r 1 , r 2 , . . . , r N ), where r i denotes the position of the ith fermion, necessarily satisfies Ψ = 0 if r i = r j for any i, j ∈ {1, 2, . . . , N }. Note that each coordinate r i corresponds to a fermion which is either spin-up or spindown, and this must be considered when calculating the explicit wavefunction. However, for the sake of simplicity, we do not incorporate explicit references to spin in the notation. Then, the wavefunction is proportional to the anti-symmetric Slater determinant, (46) where S N is the symmetric group defined on a set of size N and hence σ is any permutation of the elements {1, 2, . . . , N }, sgn(σ) is the parity of the permutation (+1 if it takes an even number of pairwise swaps to return the numbers into the original order 1, 2, . . . , N and −1 otherwise) and ϕ i (r σi ) is the quantum harmonic oscillator eigenfunction defined in Eq. (45). The general wavefunction for an N -fermion system is then expressed as [107] Ψ(r 1 , r 2 , . . . , r N (47) where we sum over all N ! permutations of the coordinates {r 1 , r 2 , . . . , r N }, and θ(r 1 , r 2 , . . . , r N ) = 1 if r 1 < r 2 < · · · < r N and 0 otherwise and a k are simple coefficients. This solution was, for example, given explicitly for the N = 2 + 2 case in Volosniev et al [107]. Note here we have the simple configuration in which the net spin is zero.
From this exact wavefunction, we can analytically calculate the density and pair-correlation function for Nfermion systems in the presence of infinitely strong repulsive interactions, thus permitting us to precisely assess the accuracy of the RDM methodology in capturing observable quantities in addition to the ground-state energy in the strongly-interacting regime. These results will be discussed in the next section, with the explicit densities and pair correlations, given in Fig. 4 and Fig. 5.

VI. RESULTS
The RDM methodology equips us with the groundstate energy as well as the linear expansion coefficients of the 2-RDM in the given basis. As such, not only can we compare the ground-state energy to other, more established methodologies in the field, but we can also construct other physical observables such as the singleparticle density and the pair-correlation function. In this section, we first plot the ground-state energies of various ultracold, one-dimensional systems and compare to other methodologies, and then we investigate the singleparticle densities and pair-correlation functions for different interaction strengths in both the attractive and repulsive regimes. Moreover, we extract the eigenvalues of the 2-RDM which permits us to analyse any pairing of fermions of opposite spin in the system. Note that in all plots presented the rank, K, refers to the spin-orbital basis rank.

VI.1. Ground-State Energies
Here, we adopt the notation of N = N ↑ + N ↓ ; thus, for example, the N = 1 + 1 system is the two fermion system with one spin-up and one spin-down fermion.
In Fig. 1, we consider the most simple system that demonstrates the efficacy of the RDM method where we plot the ground-state energy, relative to the energy of the non-interacting state E F , of the N = 1 + 1 system as a function of the interaction strength g 1D . A direct comparison between the FCI energy (solid lines) and the RDM energy with the D, Q, and G conditions enforced (circles and diamonds) is made, and it is seen that the agreement is excellent, with the magnitude of difference between the FCI and RDM energies being of the order of 10 −3 ω for all values of g 1D given, as we can see in Tab. II in Appendix A. This is expected, as the D, Q, and G conditions on the 2-RDM for an N = 2 system constitute a complete set, and the resulting energy is exact, relative to the finite basis set, with the 2-RDM completely capturing all correlation in the system.
For reference, the analytic solution for the N = 1 + 1 system, as derived by Busch et al [8], is shown by the black solid line. Naturally, we expect a difference in the predicted ground-state energy of the analytic and variational results for two reasons: firstly, we are approximating a function, i.e., the exact wavefunction, which is defined on an infinite-dimensional Hilbert space as a linear expansion of a finite number of basis functions and, secondly, the exact wavefunction possesses pathological behaviour such as cusps which are not aptly modelled the continuous basis functions. Here, a comparison to the analytic results is not the objective, but rather a comparison to the exact result in the finite basis which is the result given by the FCI calculation and the match, in this case, between the FCI and RDM energy is excellent.
In Fig. 2 we plot the ground-state energy, relative to the ground-state energy of the non-interacting system, E F , for the N = 2 + 2, 4 + 4, 5 + 5, and 6 + 6 systems. We make the comparison between the energy found by the RDM method (with the D, Q, and G initially enforced and then the T1 and T2 enforced also) and the FCI energy in the N = 2 + 2 case (panel a)). We notice that the RDM method captures the energy accurately well into the strongly-interacting regime. Moreover, for the remaining systems, we compared the RDM energies to the CCSD and CCSD(T) energies and also see that the match between the CCSD(T) and RDM energy with all constraints up to the T2 condition is excellent. Such an excellent agreement has also been demonstrated for various atoms and molecules [58]. Explicit numerical data detailing this agreement is given in Tab. II and Tab. III in Appendix A.
The systems exhibited here, particularly the N = 5 + 5 and N = 6 + 6 systems, exhibit strong correlation due to the number of interacting fermions and this is augmented as the interaction strength becomes very large. However, the RDM methodology accurately captures this correlation energy, as compared to the gold-standard CCSD(T) method. This utility of the RDM method, in addition to accurately capturing the correlation in these relatively large systems, is also manifest in the fact that the fermion number, N , enters the variational procedure merely as a parameter and has no bearing on the computational complexity, which only scales with rank. Thus, we see that the RDM scales well in its ability to capture strong correlation in larger systems with no additional computational resources required (if the calculations being compared are completed with the same spin-orbital basis rank, K).
A comparison between the ground-state results in each of these systems can be made by plotting the energy, E/E F , where E is the ground-state energy and E F is the energy of the non-interacting system, against a rescaled interaction strength, defined by γ := πg 1D / √ N . Such an approach was given by Grining et al [35] wherein it was shown that the energy difference between the N = 1 + 1 system and the rescaled energy in the thermodynamic limit of an infinite number of fermions was extremely close. In Fig. 3 we plot the rescaled energy against the rescaled interaction strength. We notice the significant overlap between the energies for each of the systems considered, which is consistent with the results given by Grining et al [35]. Such an agreement extends well into repulsive and attractive interaction regimes.
Since CCSD methodologies have been shown to accurately reproduce the correlation energies in these systems [14,35] and since the RDM method matches significantly well with the CCSD results, we can couple the ground-state energies given in the panels in Fig. 2 and the rescaled energies in Fig. 3 to decisively state that the correlation in these strongly correlated systems is accurately captured by the T1 and T2 N -representability conditions in the RDM method.

VI.2. Densities, Pair-Correlations, and Eigenvalues
In addition to the ground-state energy, we can also analyse any other physical observable as they all depend on the elements of the 1-and 2-RDMs. The density of The ground-state energies (relative, in each case, to the corresponding energy of the non-interacting system, EF ) of the a) N = 2 + 2, b) N = 4 + 4, c) N = 5 + 5, and d) N = 6 + 6 systems. We compare the energies found through the RDM method with the D, Q, and G conditions simultaneously enforced (circles) and then demonstrate the shift in energy found when the additional T1 and T2 conditions are enforced (black '+' signs). In the N = 2 + 2 system (panel a)) we utilise the FCI energy as a benchmark to compare the accuracy (and also include the energy without the T2 condition enforced) and in the remaining panels we used the CCSD and CCSD(T) energies as a comparison. Examples of the numerical differences are given in Tab. II  and Tab. III in Appendix A. the physical systems, which is defined as the diagonal elements of the 1-RDM, can be found from the elements of (1) D i j from Eq. (6). The density profiles of these ultracold systems can be measured experimentally and are key observables as they define the physical extent of the interacting system and are a key indicator of the nature and strength of the interactions in the system.
In Fig. 4 we plot the density profile of the N = 2 + 2 (panel a)) and N = 6+6 (panel b)) systems, respectively. In both, the spin-orbital basis rank is K = 40 and the D, Q, G, T1 and T2 conditions are all applied. We find the density profiles for a range of different interaction strengths, from the attractive regime at g 1D = −5 to the repulsive regime at g 1D = 5, 10, and 50. We also plot the analytic density in the limit that g 1D → ∞ which we can derive directly from the analytic wavefunction in Eq. (47) by integrating over the modulus squared of the wavefunction where the proportionality accounts for the normalisation. The dotted line indicating the analytic density in Fig. 4 matches well with the density found by the RDM method at g 1D = 50, indicating that such a physical observable is aptly described by the RDM method for both the N = 2 + 2 and N = 6 + 6 systems, even in the relatively small basis rank of K = 40.
The range of interaction strengths assessed in these density plots is sufficient to observe the impacts of attractive or repulsive interactions have on the system. We see a thickening of the tails as g 1D attains larger, positive values, which is to be expected in a system with strong repulsive interactions. This is contrasted with the density profile in the attractive interaction regime, with the width being smaller and the system attaining a high peak density closer to the centre of the trap. Such effects are clearly more pronounced in the N = 6+6 system relative to the N = 2 + 2 system.
Another key physical observable is the pair-correlation function, which encodes all of the information about the pair-wise interactions the system exhibits. We define this generally as the diagonal element of the 2-RDM, There is inherent ambiguity in the definition of (2) ρ in this case, as the 2-RDM (2) D(x 1 , x 2 ; x 1 , x 2 ) generally encodes information about the correlation between particles of the same species (which do not occur in this system due to Pauli exclusion and the interactions being pointlike) as well as between different species. Since we are considering contact interactions, it is the latter components of the 2-RDM which we consider and thus we only sum over the appropriate indices in Eq. (7) to garner the pair-correlation function. As in the case of the densities, we can also calculate the analytic pair-correlation function in the limit g 1D → ∞ directly from the definition of the wavefunction in Eq. (47). In Fig. 5 we plot the 'anti-diagonal' component of the pair-correlation function, (2) ρ(x, −x), for the N = 2 + 2 (panel a)) and N = 4 + 4 (panel b)) systems, for a range x (a ho ) 4. The density, which is equal to the diagonal element of the 1-RDM, ρ(x) = (1) D(x; x), for the N = 2+2 (panel a)) and the N = 6 + 6 (panel b)) systems, plotted as a function of x for four different interaction strengths, g = −5, 5, 10, and 50. In all cases K = 40 and the D, Q, G, T1 and T2 conditions are all applied. The density is normalised to the fermion number N . The density derived from the analytic wavefunction in the presence of infinitely strong interactions is given by the black dotted line in both cases.
of values of the contact interaction strength, from the attractive regime at g = −5 into the repulsive regime at g = 5, 10, and g = 50. We also include the analytic pair-correlation in the limit that g 1D → ∞, which is indicated by the dotted black line in both cases. We notice the ability of the RDM methodology in capturing the pair-correlation in the strongly interacting limit by comparing the RDM pair-correlation at g 1D = 50 to the exact solution at g 1D → ∞.
In all cases, the spin-orbital basis rank is K = 40 and the D, Q, G, T1 and T2 conditions are simultaneously enforced. The pair-correlation function, as defined as x (a ho ) 5. The pair-correlation function, (2) ρ(x, −x), for the N = 2+2 (panel a)) and N = 4+4 (panel b)) systems, plotted as a function of x for four different interaction strengths, g = −5, 5, 10, and 50. In all cases r = 40 and the D, Q, G, T1 and T2 conditions are all applied.
the diagonal element of the 2-RDM, is a probability distribution function for two particles (of opposite spin, in this case) occupying positions r 1 and r 2 , and hence the anti-diagonal component thereof is readily interpreted as a probability distribution function for two opposite spin particles to be occupying positions equidistant from the centre of the trap. As we can see, the attractive interactions significantly increase the likelihood with which two particles will be found in mutual proximity (near the centre of the trap).
A final characteristic we consider in this work is the spectral decomposition of the 2-RDM. The eigenvalues of the components of the 2-RDM corresponding to the interaction between particles of opposite spin gives us insight into the existence of pairing in the system. The archetypal example of pairing in fermionic system is the existence of Cooper pairs in low-temperature electronic systems [110], which is the driving mechanism for the superconductivity as explained in the Bardeen-Cooper-Schreiffer (BCS) theory [111,112]. However, other pairing mechanisms in ultracold system have been studied both theoretically and experimentally [113][114][115][116][117]. The ultracold systems analysed in this work are known to exhibit pairing [118] and here we demonstrate that the RDM methodology, in addition to capturing the salient features of these systems as already discussed, accurately predicts the occurrence of such a phenomenon.
In Fig. 6 we plot the first eight eigenvalues of the 2-RDM as a function of the interaction strength in the attractive regime for both the N = 5 + 5 (Fig. 6 a)) and N = 6 + 6 ( Fig. 6 b)) systems. In both cases, the first-order eigenvalue for each interaction strength is indicated by the circles and the remaining seven are indicated by the squares. We note that a single eigenvalue dominates above the others, indicating a BCS-like pairing manifesting for two opposite-spin particles occupying the same spatial orbital [119]. Such results are in agreement with previous results demonstrated in the literature [118] and, moreover, are consistent with the behaviour of the pair-correlation function for attractive interactions as in Fig. 5, in which the likelihood of two opposite-spin particles being proximate to one another is significant.

VI.3. A Note on Efficiency
A key facet of the effectiveness of the RDM method is the relative lack of memory and time requirements when it comes to computational implementation compared to the FCI method, as the RDM method scales polynomially in the basis rank and does not scale with particle number at all (with the FCI method scaling exponentially in terms of the basis rank and the particle number, at best [120]). The computational complexity of the interiorpoint methods scaled approximately as K 16 , where K is the basis rank. Such a poor scaling meant that the method was confined to small atomic and molecular systems. However, faster algorithms, such as those used by SDPNAL+ result in a scaling of approximately K 6 to K 8 . Such an improvement in the scaling of the complexity was also manifest in the boundary-point method developed by Mazziotti [54], which demonstrated a tento twenty-fold increase in efficiency to other calculations at the time.
To demonstrate the efficiency, we found, in Tab. I, the CPU time (in minutes) for various calculations, and the dependency of the CPU time on fermion number and basis rank for the 2-RDM and FCI calculations, as well as the percentage of the correlation energy captured by the 2-RDM method, with the FCI energy being defined as 100% correlated. The 2-RDM method does not scale with fermion number, which is an enormous advantage over methods that require to calculation of the N -fermion wavefunction. We can see the large increase in CPU time for large rank and fermion number for the FCI calculations, with the relatively slow increase with rank of the  Enforcing further constraints on the 2-RDM, such as the T1 and T2 constraints results in a much larger SDP and hence a larger CPU time. Generally speaking, the largest hindrance to the widespread adoption of this method was the poor scaling of the computational complexity with basis rank as one increased the amount of N -representability conditions enforced. However, recent effort has significantly reduced this scaling from approximately K 9 to K 6 [121,122]. I. CPU time (in minutes) for various systems and spin-orbital basis ranks for the 2-RDM and FCI calculations. All 2-RDM calculations simultaneously enforced the D, Q and G conditions. We note that the 2-RDM calculations scale only as a function of the rank K, and hence only one calculation need be performed for each, whereas the FCI method scales exponentially with both basis size and rank, as seen for each rank as we increase the particle number. Calculations were performed on a laptop with Intel Core i7 2.60GHz processors (64GB RAM) using 6 cores. The accuracy of the sample 2-RDM calculations is indicated by the percentage of the correlation energy (where the correlation energy is defined by Ecorr = E HF − E FCI , i.e., the difference between the mean-field HF energy and the exact FCI energy). In all calculations, g 1D = 10.

VII. CONCLUDING REMARKS
In this work, we outlined the theoretical development of the RDM method, as well as briefly demonstrating how to complete a simple implementation via a formulation as an SDP. Such a method has seen wide success in atomic and molecular applications, and here we extend the range of applicability by considering a simple, highly correlated few-fermion system. Specifically, we applied the RDM method with the D, Q, G, T1, and T2 N -representability constraints to a trapped, onedimensional system of fermions interacting via a contact interaction, while constraining the system to have an equal number of spin-up and spin-down fermions. The ground-state energy was found for N = 2, 4, 8, 10 and 12 particle systems and was found to match well, for a range of interaction strengths in the positive and negative regimes, with the best theoretical results in the given harmonic oscillator, basis, provided by the FCI or exact energy for the N = 2 and N = 4 system, and by the CCSD and CCSD(T) method for the N = 8, 10, and 12 system.
Specifically, for the N = 1 + 1 case, the ground-state energy given by the RDM method with the DQG conditions implemented agreed with the exact FCI energy for all interaction strengths, −5 ≤ g ≤ 20, which is expected. In the N = 2 + 2, 4 + 4, 5 + 5 and N = 6 + 6 cases, it was demonstrated that while constraining the 2-RDM with the D, Q, G, T1, and T2 conditions, strong correlation was accurately captured, with the ground-state energy being close to that given by FCI and CCSD(T) methods. This demonstrates the effectiveness of the RDM method when it comes to capturing the ground-state energies in strongly correlated systems, and further demonstrates that this method is an extremely promising alternative to typical methods based on the complete knowledge of the N -fermion wavefunction. Moreover, we demonstrated that the three-body N -representability constraints are sufficient in a variational procedure to minimise the 2-RDM, and higher-order constraints need not be considered.
The density and components of the pair-correlation function we also considered for the N = 2 + 2 and N = 6 + 6, and N = 2 + 2 and N = 4 + 4 systems, respectively, with a basis rank of K = 40 and simultaneous enforcement of the D, Q, G, T1 and T2 conditions. These quantities are calculated directly from the variationally determined components of the 1-and 2-RDMs. Such observables were found in the strongly interacting regime, where g 1D = 50, and were found to match well with the anayltic result for the observables in the limit of infinitely strong repulsive interactions. By means of Eq. (6) and Eq. (7) any physical observables can be readily expressed in terms of the components of (1) D and (2) D. We also demonstrated that important features of these systems such as the phenomenon of BCS-like pairing was manifested in the spectral decomposition of the 2-RDM. As such, the key features of these systems are all found naturally by considering the variationally determined elements of the 1-and 2-RDMs.
The strong correlation in these one-dimensional systems, which occurs in the regime of |g 1D | 1, was accurately captured by the RDM method after implementing the three-body constraints known as the T1 and T2 conditions. This shows that, in addition to the RDM method accurately capturing strong electron correlation in atomic and molecular systems, the RDM method accurately accounts for strong correlation in ultracold fewfermion systems where g 1D can be made arbitrarily high. An important point to note, with regard to this method, is that the introduction of higher-order constraints provides a mechanism for understanding the limitations of lower order calculations. Additionally, for a given rank, utilisation of the RDM method in parallel with conventional wavefunction techniques can allow one to evaluate an upper and lower bound for the ground-state energy of the system.
The choice of basis, namely the harmonic oscillator basis, naturally hinders the attempt to find extremely accurate ground-state energies relative to the known, analytic case for the N = 1 + 1 system, and other approach for larger N . However, what has been clearly demonstrated is that the RDM method gives accurate results relative to the FCI and CCSD(T) methods for relatively large systems and strong interaction strengths in the same basis.
The accuracy of the RDM method for these simple systems opens the door to applying the same method for more complex interacting systems. An extension to the analogous system studied here in two and three di-mensions is a natural case for further investigation. Further extensions using this methodology could include spin imbalanced systems, where N ↑ = N ↓ , as well as massimbalanced systems. Minor imbalances in these systems can be exaggerated to permit the study of single impurities in otherwise homogeneous systems [123][124][125]. Alternate interactions, such as soft-core or dipolar interactions can be analysed with this method also. Such changes amount to mere substitutions in the Hamiltonian matrix elements, and modifications of the particle number and net spin, which are simple parameters in the calculation and requires no modification of the methodology as a whole nor any increase in computational complexity. This demonstrates the universality of the RDM approach to many-body quantum systems with pairwise interactions and exemplifies why it is a promising method in the field.
The 2-RDM approach adopted here achieves the goals of N -representability and universality that have driven the development of DFT, which is the 1-RDM method for fermion systems. The exponential complexity of Nfermion wavefunction-based approaches has been circumvented in favour of a method with polynomial scaling characteristics. The systematic order-by-order application of N -representability constraints leads to a computational scheme for 2-RDM that can be made as accurate as required; it is a "gold standard" comparable to full configuration interaction (FCI). Nevertheless, the computational cost of 2-RDM exhibits an intrinsic computational complexity far greater than that of DFT, and comparable with coupled-pair methods, such as the CC approximation. The 2-RDM method is, however, more than just a convenient computational tool. It provides direct insights into the relationship between the energy of a fermionic system and its 2-RDM through the simple process of integration over the coordinates of a single fermion. The physical principles that underpin the properties of N -fermion wavefunctions that are used to derive the Hohenberg-Kohn theorems of DFT apply with equal validity to an N -representable 2-RDM. This suggests the potential for a parallel development of 1-RDM methods based on the availability of 2-RDM approaches for complex fermion systems.  TABLE II. Comparison between the energy found from the FCI calculation, E FCI , and the energy found by the RDM method, with the RDM energy given as a percentage of E FCI . Here, we consider the N = 1 + 1 and N = 2 + 2 systems for g 1D = −5, 10, and 20.
energy (as % of E CCSD(T) ) g 1D N E CCSD(T) E DQG E DQGT1 E DQGT1T2 4 + 4 -5.61208 95 Comparison between the energy found from the CCSD(T) calculation, E CCSD(T) , and the energy found by the RDM method, with the RDM energy given as a percentage of E CCSD(T) . Here, we consider the N = 1 + 1 and N = 2 + 2 systems for g 1D = −5, 10, and 20.
energy captured by the RDM method as a percentage of the energy captured by the FCI (for the N = 1 + 1 and N = 2 + 2 systems) and the CCSD and CCSD(T) methods (N = 4 + 4, 5 + 5, and 6 + 6) for the interactions of strength g 1D = −5, 10, and 20.
The ground-state energy found by the RDM method is lower, in all cases, compared to the FCI or CCSD/CCSD(T) energies. This is indicative of the fact that if one adds successive N -representability constraints to the variational determination of the 2-RDM we approach the ground-state energy from below. Thus, lower values of the percentage of the energy captured by the RDM method (when compared to the FCI or CCSD/CCSD(T) which we take as a benchmark) indicate that further, higher-order constraints are needed. It must be noted, however, that the numerical results further demonstrate the ability of the RDM method to capture the strong correlation, even in the relatively large N = 5 + 5 and N = 6 + 6 systems.