Variational ansatz for $p$-wave fermions confined in a one-dimensional harmonic trap

We propose a very accurate and efficient variational scheme for the ground state of the system of $p$-wave attractively interacting fermions confined in a one-dimensional harmonic trap. By the construction, the method takes the non-analytical part of interactions exactly into account and thus it approximates the true ground-state wave function in a whole range of interactions very accurately. Within the method, we determine different properties of the system for a different number of particles and different interactions. In this way, we explore how the system and its features transit from the ideal non-interacting Fermi gas to the system of infinitely strong attractions. The presented method of including zero-range interactions is very universal and may be easily generalized to other one-dimensional confinements.


Introduction
Appropriate description of strongly correlated quantum many-body systems offering an adequate explanation of their different measurable properties is one of the most challenging tasks for theoretical physics from over 60 years [1]. Besides theoretical reasons, it became fundamentally important recently due to tremendous progress in quantum engineering giving opportunities to coherent control of matter and light on atomic scales where an accurate theoretical description is required. The fundamental obstacle for all the straightforward descriptions of many-body systems originates in the mathematical complexity of a many-body Schrödinger equation which, in fact, can be analytically solved only in several specific cases. To the most famous examples belong: the Moshinsky model and its variations [2,3,4,5,6], the Lieb-Liniger model [7,8], or the Calogero-Sutherland model [9,10]. Even in the case of only two interacting particles the list is not significantly extended and contains only few additional solutions: the famous Busch et al. solution for s-wave contact forces [11], its variations [12,13,14] and generalization to p-wave forces [15], the Gao solutions for 1/r 6 and 1/r 3 potentials [16,17], and specific solutions for finite-range interactions modeled by a step function [18,19].
The situation becomes even more complicated and almost hopeless when interparticle interactions are modeled by non-smooth functions. Then, many direct numerical attempts are not able to capture subtle features of such interactions and simply break down. One of the simplest, but still realistic models dominated by these kinds of problems is the one-dimensional model of N identical spinless fermions interacting via two-body zero-range forces [20,21]. In contrast to bosonic systems, in this case, the Pauli exclusion principle precludes any scattering in the s-wave channel and the first non-vanishing contribution to interactions comes from the p-wave scattering [22,23]. In one-dimensional case this interaction acts highly counterintuitively: if φ 1 (x) and φ 2 (x) are two different wave functions describing relative motion of two fermions then the matrix element of p-wave interaction is proportional to the product of their spatial derivatives at the origin, Bypassing this difficulty in any numerical treatment is not an easy task if the mentioned non-analytically of many-body wave functions are not taken exactly into account. Importantly, special care on this problem needs to be put when then the number of particles is not large since then any inaccurate approximation may lead to significant discrepancies in physical predictions. In this work, based on our previous experience with zero-range forces in the bosonic case [24], we propose a very accurate and very efficient way to find approximate ground-state wave function of a few p-wave interacting fermions. To show that the method proposed is indeed profitable, we focus on the generic problem of N identical fermions attractively interacting via p-wave forces confined in a one-dimensional harmonic trap. Preliminarily, the model has been studied already, mostly in the limit of infinite attractions [25,26,27,28]. We show that within our framework we can determine with high accuracy not only the groundstate energy and single-particle properties of the system (determined previously for other scenarios [29,30,31]) but also we can go much further and study the interparticle correlations in position and momentum domain. As an example, we display two-particle probability densities highlighting intriguing correlations, especially in the momentum domain. Importantly, the variational method proposed opens a route to capture different properties of the system for intermediate interactions and therefore to observe their evolution when interactions are tuned along a whole range.

The system
In our work we consider the system of N identical spinless fermions of mass m confined in a one-dimensional harmonic trap of frequency ω and interacting with zero-range interactions. For convenience, in the following considerations we will express all quantities in natural units of the problem studied, i.e., all energies, lengths, and momenta will be expressed in units of ω, /mω, and √ mω, respectively. As already mentioned in the introduction, the s-wave scattering between particles is not present and the first non-vanishing contribution to interactions comes from the zerorange p-wave interactions. In a one-dimensional geometry, the interaction potential can be expressed formally as [32,23,20,33] where x is a relative distance between particles and g F is the effective p-wave interaction strength. Thus, the Hamiltonian of the system studied has a form It turns out that in the one-dimensional scenario the fermionic Hamiltonian (2) is exactly equivalent to the problem described by the non-interacting Hamiltonian provided that the many-body eigenstate wave function ψ(x 1 , . . . , x N ) is antisymmetric under an exchange of any two positions and it additionally supports the contact condition for each pair of fermions of the form [21,34,30] ψ( where x ij = x i − x j . This equivalence simply means that a whole effect of p-wave interactions between particles is encoded directly in the condition (4). Our aim is to find convenient approximate form of the ground-state wave function fulfilling this condition exactly.

Variational approach
Before we present our construction of the trial function for the N particle system, let us note that the problem studied has known analytical solutions in the case of N = 2 fermions [15]. The eigenenergies and corresponding eigenvectors are found analogously as in the celebrated Busch et. al problem of two s-wave interacting boson [11]. For p-wave fermions all the wave functions are expressed in terms of the confluent hyperbolic function U (a, b, z) as where the parameter is determined by the contact condition (4). For fixed interaction strength g F this condition means that the wave functions (5) are egienstates of the Hamiltonian (2) provided that corresponding parameters are successive solutions of the following transcendental equation: Eigenenergies are directly related to via the relation E n = + 2n + 1. Consequently, the ground-state solution corresponds to the state with the lowest . The exact solution of the problem with more than two interacting fermions is known only in the two extreme limits, i.e., the non-interacting case (g F = 0) and infinitely strong attractions (g F = −∞). It should be noted however, that for any number of particles N and for any attractive interaction strength g F there exists a direct and rigorous mapping of the p-wave fermions ground-state wave function ψ(x 1 , . . . , x N ) to the ground-state wave function of s-wave repulsively interacting bosons Φ(x 1 , . . . , x N ) described by the many-body Hamiltonian of the form [35,36]: where g B is the effective s-wave interaction strength between bosons. Let us formulate the mapping more precisely. If Φ(x 1 , . . . , x N ) is the ground-state wave function of the Hamiltonian (7), then the wave function defined as is the exact ground-state wave function of the Hamiltonian (2) describing p-wave fermions provided that the corresponding interaction strengths g B and g F fulfill the condition g B · g F = −2. Consequently, from the principal point of view, the question of finding the ground state of p-wave fermions addressed here is rigorously equivalent to another problem of interacting bosons. However, the problem of efficient construction of the wave function remains unsolved since the exact form of bosonic wave function is not known. Therefore, is it still interesting to find accurate and efficient method of constructing many-body wave function for interacting p-wave fermions. One of the natural method of finding the many-body ground-state wave function is to propose a reasonable family of variational trial functions appropriately tailored to capture the most important features of the system. Having in hand the exact twobody solution (5) which is valid for any interaction strength, one can introduce the Jastrow-like trial function for the problem of N particles [37]. In this approximation it is assumed that the most prominent part of inter-particle correlations are captured in the two-body sector described by the two-body solution. Therefore, the trial wave function has a form: with a rescaled two-body solution In this approach, α is the variational parameter of the family and α are determined by solving the equation In principle, by minimizing the energy functional one can obtain approximation to the ground-state of interacting system of N fermions and corresponding energy. However, in practice this approach is very demanding due to a quite large computational complexity. This complexity originates in a tangled definition of the rescaled two-body solution ψ α (x 1 , x 2 ) through the hypergeometric function determined by α being a solution of transcendental equation (6). Consequently, the method cannot be easily used for large number of particles. A very similar problem was revealed recently in the case of interacting bosons [24]. To overcome the difficulty described above, instead of Υ α (x 1 , . . . , x N ), we propose to use different variational ansatz which turns out to be reasonably accurate and numerically very efficient. This approach is still based on the assumption that the dominant part of correlations has two-body origins and can be written as However, the correlated pair function φ α (x) is significantly simplified to the form where α plays a role of the variational parameter. Importantly, the two-body trial function φ α (x) has appropriate properties in the vicinity of x = 0: Thus, the whole variational wave function (13) fulfills automatically the discontinuity condition (4) and therefore it includes exactly all effects of inter-particle interactions. Moreover, the ansatz reproduces rigorously the system's ground state in the two extreme limits for any number of particles N . First, in the limit of non-interacting fermions (g F → 0) the parameter α → 0 and as a result

Ground-state energy
The variational trial wave function (13) is a very efficient tool for determining different properties of many interacting p-wave fermions for intermediate interaction strengths. First, it allows us to find quickly the ground-state energy for a given number of particles N . As shown in Fig. 1, ground-state energies obtained in this way nicely interpolate between two extreme points -the non-interacting energy E 0 (N ) = N 2 /2 and the ground-state energy in the limit of infinite attractions E −∞ (N ) = N/2. Moreover, in the case N = 2, the energies almost ideally match the exact values obtained from relation (6) (red dots). To show that the method can be successfully used for a quite large number of particles, in the right panel of Fig. 1 we plot the ground-state energy as a function of the number of particles for some particular values of interactions.

One-body properties
When the variational ground-state energy is found one automatically has an approximate representation of the many-body ground state for interacting p-wave fermions. This enables one to study different properties of the system. All singleparticle ones are encoded in the single-particle reduced density matrix. In the position representation it can be calculated straightforwardly from the ground-state wave function as: Its diagonal part n (1) (x) = ρ (1) (x, x) encodes the single-particle density profile. Taking the case of N = 4 particles as an instructive example, in the first two rows in Fig. 2 we plot these quantities for different interactions. It is clearly seen that along with increasing p-wave attractions, the density profile and single-particle density matrix change their shapes. In the limit of infinite attractions (g F = −∞) both of them resemble corresponding features of non-interacting bosons which is in a full accordance with the mapping mentioned above. The situation changes significantly when, instead of spatial, momentum properties of the system are discussed, since momentum distributions cannot be easily Single-particle properties of the system of N = 4 fermions and different interaction strengths. First two rows correspond to the single-particle density profile n (1) (x) and the single-particle reduced density matrix ρ (1) (x, x ) in the position representation, while next two rows to the corresponding quantitiesñ (1) (p) and ρ (1) (p, p ) in the momentum representation. Corresponding results for N = 3 and N = 5 are displayed in Fig. A1 and Fig. A2, respectively. determined by simple mapping from the associated s-wave bosonic system. This unfeasibility originates in a simple fact that the mapping between attractive fermions and repulsive bosons is a non-unitary transformation in the position domain.
The simplest properties of the system in the momentum domain are encoded in the single-particle density matrix defined as and its diagonal part (the momentum density profile) of the form In two bottom rows in Fig. 2 we display these distributions for N = 4 and corresponding interactions. It turns out that along with increasing attractions in the system higher momenta of single fermions are accessible. Thank to the method used, now one has an access to the ground-state distributions also for intermediate interactions.
In this way one can easily observe how the two exterior peaks present in the distribution of non-interacting system are smeared along with increasing attractions while remaining two are enhanced. From this point of view, it is also very instructive to compare this behavior for different number of particles. Therefore, in the Appendix we present them for another cases with N = 3 and N = 5 particles (see Fig. A1 and Fig. A2, respectively). As it is seen, in these cases (note odd number of particles) the central peak is enhanced while only two the most external are smeared. At this point let us also mention that easy access to the full single-particle density matrix of the system provided by the proposed variational scheme gives us also a direct way to quantify non-classical inter-particle correlations in the system. Most simply, this can be done by performing spectral decomposition of the single-particle density matrix where λ i and functions η i (x) are eigenvalues and natural orbitals of the single-particle density matrix. Since the system contains N indistinguishable fermions, even in the non-interacting case (g F = 0) the decomposition is not trivial and has exactly N non-zero eigenvalues λ 1 = λ 2 = . . . = λ N . They correspond to N single-particle orbitals forming the Slater determinant describing the state of the system. For nonvanishing interactions more than N orbitals contribute to the density matrix and the situation becomes more complicated. Nonetheless, still a general structure of the single-particle density matrix is rigorously known [38,39]. For the even number of particles N all the eigenvalues {λ i } are evenly degenerated. In contrast, for odd N , exactly one of the eigenvalues is always equal to 1/N and remaining ones are evenly degenerated. As an example, in the top right panel in Fig. 3 we display the dependence of a few the largest eigenvalues λ as functions of interactions for the system of N = 4 fermions. Corresponding plots for N = 3 and N = 5 particles are supplemented in the bottom of the same figure. It is clear, that our variational approach appropriately reproduces the structure of the reduced density matrix. Interestingly, close to the non-interacting limit, all the eigenvalues rapidly change their values and for not very strong interaction strengths (g F ≈ −4 for N = 4) saturate at values being very close to their values for infinite attractions. It may suggest that many important one-body features of the system achieved in the limit g F → −∞ are exhibited by the system already for relatively weak interactions. Obviously, the exact number of contributing orbitals and corresponding eigenvalues is not determined by general theorems and they depend on the amount of non-trivial correlations present in the system. They can be quantified by the von Neumann entanglement entropy defined as The entropy is obviously bounded from below by its value in the non-interacting limit, S 0 = ln N . In the right panel in Fig. 3 we present the dependence of the entanglement entropy S on interactions for a different number of particles. It is clear, that entropy monotonically increases with interaction strength which signals a monotonic increase of correlations in the system. However, exactly as anticipated by the behavior of eigenvalues, the entropy quickly saturates on its value reached in the limit of infinite attractions. Additionally, in the case of N = 2 fermions, we mark the exact values of the entropy provided by the exact solution (5). With this comparison, it is clear, that the variational approach proposed appropriately captures also quantitative predictions for non-trivial one-body coherence.

Two-body correlations
The most important advantage of our variational method is its ability to predict higherorder correlations between interacting fermions. Since the full many-body wave function is appropriately represented in a whole range of interactions, it gives a direct view on changes of different non-trivial correlations under tuning interaction strength. For example, one can easily consider two-particle density profiles in the position and momentum domains. They are defined as n (2) (p 1 , p 2 ) = 1 4π 2 dxdx dydy ρ (2) (x, y, x , y )e −ip 1 (x−x ) e −ip 2 (y−y ) , (21b) where we introduced the two-particle reduced density matrix of the form Physically, the profiles (21a) and (21b) can be interpreted directly as probabilities density of finding two fermions with positions (x 1 , x 2 ) or momenta (p 1 , p 2 ) in a simultaneous measurement of two particles. Therefore, they are the simplest quantities capturing geometrical (in positions or in momenta) arrangement of particles in the many-body ground state. It turns out that properties of these two distributions crucially depend on interactions. To visualize this, in Fig. 4 we plot them for the system of N = 4 particles exactly for the same interaction strengths as in Fig. 2. Interestingly, the the two-particle distribution in the position domain undergoes a smooth transition from the square-like to the circle-like shape when attractive forces are enhanced. This transition is assisted by a significant reduction of the forbidden For comparison, the exact results obtained for N = 2 case are marked by red dots. In the two limiting cases g F = 0 and g F → −∞, the mean squared distance can be determined analytically and it is equal to N + 1 and 1, respectively. Note that for better visibility, on the left panel we use a nonlinear scale for interactions. region along the diagonal x 1 = x 2 . This behavior does not qualitatively depend on the number of particles (see Fig. A1 and Fig. A2 in Appendix). Contrary, in the momentum space the situation is completely different and it strongly depends on the parity of the number of fermions. While in the non-interacting case, due to the symmetry of canonical variables x ↔ p, the correlation function has also a square-like shape, along with increasing attractions some enhancements of two-body correlations in particular directions appear. Although in the case of N = 4 particles almost only the ordinary pairing of opposite momenta is supported (clear enhancement along the line p 1 +p 2 = 0), for odd number of particles (see Fig. A1 and Fig. A2) we notice additional strong enhancement of pairs in which only one of particles carries all momentum (enhancement along lines p 1 = 0 and p 2 = 0, respectively). The dependence of the two-particle correlations on interactions and the number of particles can be also visualized when the mean squared distance between two fermions is considered. It is defined directly as The quantity reflects the spreading of the two-particle density profile in the position domain and its value is known for any number of particles N in two extreme limits. It is equal to N + 1 for the non-interacting system (g F = 0) while for infinite attractions (g F → −∞) it is independent on N and equal exactly to 1. For intermediate interactions, the mean squared distance σ 2 can be quite easily determined by the variational method proposed. In Fig. 5 we display the results obtained for a different number of particles up to N = 10 and a whole range of attractions. It is clear, that for all number of particles considered, the squared distance σ 2 rapidly decreases with attractions and it quickly achieves its asymptotic value 1. Even for a quite large number of particles, N = 10, the squared distance σ 2 is less than 2 if interactions are not weaker than g F = −4.0. This observation supports our previous singleparticle conclusions that different properties of the system with infinite attractions are revealed already for relatively weak interactions.

Conclusions
We have shown that the many-body ground-state of interacting p-wave fermions confined in a one-dimensional harmonic trap can be well-approximated by a simple variational wave function of the Jastrow type. However, in contrast to the original idea of Jastrow, instead of utilizing a known analytical solution of the corresponding twobody problem, we propose to use a much simpler correlated pair function which takes inter-particle interactions precisely into account but significantly simplifies numerical calculations. In this way, we were able to determine different single-and twoparticle properties (in the position as well as in the momentum domain) of the system containing up to 10 particles in a whole range of attractive interactions. It is worth pointing out that the scheme proposed is very flexible and general since the contact condition (4) is included in the trial wave function independently on the external confinement. Therefore, by simple modifications of the analytical part in (13) one can repeat the scheme for other one-dimensional traps. Moreover, whenever the accuracy of the approximation is insufficient, the same modification can be exploited to propose another, more precise, family of variational trial functions without modifying the pair-correlation part. Finally, since a whole strategy relies on an appropriate inclusion of the contact condition, a very similar scheme can be used for zero-range forces other than p-wave. For completeness of the discussion in this Appendix we present results obtained for different number of particles. In Fig. A1 and Fig. A2 we display the same quantities as shown in Fig. 2 and Fig. 4 but for N = 3 and N = 5 particles, respectively.