{'fraction_non_alphanumeric': 0.03018853225303833, 'fraction_numerical': 0.006427235899033967, 'mean_word_length': 4.837426102773988, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Four fermions in a one-dimensional harmonic trap: Accuracy of a variational-ansatz approach

Detailed analysis of the system of four interacting ultra-cold fermions confined in a one-dimensional harmonic trap is performed. The analysis is done in the framework of a simple variational ansatz for the many-body ground state and its predictions are confronted with the results of numerically exact diagonalization of the many-body Hamiltonian. Short discussion on the role of the quantum statistics, i.e. Bose-Bose and Bose-Fermi mixtures is also presented. It is concluded that the variational ansatz, although seemed to be oversimplified, gives surprisingly good predictions of many different quantities for mixtures of equal as well as different mass systems. The result may have some experimental importance since it gives quite simple and validated method for describing experimental outputs.

The physics of few quantum particles is extremely difficult to be analyzed without any approximations. It comes from the simple observation that 'a few' is too many to use a straightforward method for one-and twobody physics, and at the same time it is still not enough to adopt methods of statistical many-body theory and mean-field description [35][36][37]. Therefore, one has to find completely different approaches to the problem (for example those, which were up to now in a domain of nuclear physics [38]). Independently of these facts, there always exists a temptation to describe complicated few-body problem with evidently oversimplified methods. One of these kind of approaches is based on different implementations of the variational-ansatz method.
In this paper we want to investigate the properties of a system of four fermionic atoms confined in a onedimensional harmonic trap obtained via a simple variational method and validate these results. The method is based on an assumption that the ground state of a many-body interacting system can be almost perfectly superposed from two limiting many-body states, i.e., the ground states obtained for vanishing and very strong repulsions [39]. Since the method was successfully adopted for systems of two and three quantum particles (and for a particular class of polaron systems with up to six bodies), a natural question about the validity of this assumption for larger number of particles arises with other system compositions. Here we try to answer this question by comparing predictions of the ansatz with predictions of numerically exact diagonalization of the fourbody Hamiltonian. A comparison is done on various levels by considering many different quantities that, in principle, may be extracted from the experimental data. We stress that we consider the experimentally relevant situation where the particles in our system have different masses. This is a particularly difficult issue for onedimensional systems. Such systems cannot be addressed using for instance the Bethe ansatz as mass differences will generically break the assumption of non-diffractive scattering. This assumption is central to the traditional Bethe ansatz approach to generate exact solutions of onedimensional many-body systems [2,3]. This implies that a simple approach to mass imbalanced systems is highly desirable.
In Sec. II a brief description of the system under study is given and both complementary methods of treatment, i.e., the interpolatory ansatz and the exact diagonalization are briefly characterized. In Sec. III we compare different predictions of both methods and we discuss disclosed discrepancies. Finally, in Sec. IV we give some remarks on four-body systems with other quantum statistics, we discuss some possible extensions of the variational method and conclude briefly.

II. THE MODEL
The system studied.-We consider N a = N b = 2 fermionic particles confined in an external onedimensional harmonic potential of frequency ω. In principle, the particles of different kinds may have different masses, i.e., m a = m b . We assume that interactions between particles can be described with two-body contact δ-like potential. In this case, due to the fermionic nature of particles, the interactions are present only between particles of different components. The Hamiltonian of the system has the form Notice that since the masses are allowed to be different, the Bose-Fermi mapping [40][41][42] cannot be applied. However, in extreme limits the exact ground-state wave function is always known. For vanishing interactions, g = 0, particles occupy the two lowest single-particle orbitals of corresponding harmonic oscillators. In the case of infinitely strong interactions, 1/g = 0, a semi-analytical expression for the exact four-body ground state was found recently in [43] using the methods introduced in [44,45]. In general, the properties of the ground state for intermediate interactions cannot be found analytically and one needs to use numerical or approximate methods. Interpolatory ansatz.-Quite recently, it was proposed to use a very simple variational method based on the assumption that the ground state of the system for any interaction can be well approximated by an appropriate superposition of the ground states in the limiting cases: The coefficients α(g) and β(g) are determined by minimizing an expectation value of the many-body Hamiltonian (1) in this state. Note, that the many-body states |Ψ 0 and |Ψ ∞ are not necessarily orthogonal. Therefore, the variational parameters fulfill non-natural normalization conditions. The detailed prescription for obtaining appropriate variational parameters α(g) and β(g) was discussed in [39]. For the completeness of our discussion we include a brief discussion of the method in the Appendix. A small modification of the method also mentioned in the Appendix which substantially improves predictions of the ground-state energy is discussed in further analysis.
Although the ansatz seems to be highly oversimplified it was used for systems with equal masses with surprisingly good results. Here we want to make a comprehensive study of the accuracy of the ansatz when different quantities and interparticle correlations extracted from the ground state are considered. Especially, we are interested in the cases when the particles belonging to the different components have different masses. To find quantitative answers to this open questions we perform the numerically exact diagonalization of the many-body Hamiltonian Eq. (1), we find its exact ground state as a function of interactions and we compare different quantities with predictions of the variational ansatz.
Numerical diagonalization.-The exact diagonalization is performed in a straightforward and well-established way. First, we express the many-body Hamiltonian Eq. (1) in a matrix form in an appropriate Fock basis. It can be done by expressing all many-body states of the system in the basis composed as products of singleparticle orbitals: where ϕ a,k (q) are eigenstates of corresponding singleparticle harmonic oscillators, i.e., and A . is the anti-symmetrization operator in the appropriate subspace of indistinguishable fermions assuring that: |kl; mn = −|lk; mn = −|kl; nm .
Assuming some sufficiently large cutoff k ≤ N max of the considered single-particle excitations one can calculate all matrix elements of the Hamiltonian Eq. (1). The resulting matrix is diagonalized to find the exact ground state of the system |Φ(g) and its energy E(g). In our case, the exact diagonalization is performed with the Arnoldi method [46] that was used previously with a great success for similar models [47][48][49][50][51]. Alternative diagonalization routines that exploit effective interactions are also very efficient for all interaction strengths [52,53], although these methods have yet to be extended to the case with particles of different mass. In the following, many-body wave functions in position representation corresponding to states |Ψ(g) and |Φ(g) will be denoted as Ψ g (q 1 , q 2 ; q 3 , q 4 ) and Φ g (q 1 , q 2 ; q 3 , q 4 ), respectively. Additionally, we introduce a dimensionless parameter µ = m a /m b for the mass ratio of atoms from different components.

III. QUALITY OF THE ANSATZ WAVE FUNCTION
The ground-state energy.-The quality of the assumed form of the variational wave function can be examined in various ways depending on the physical quantity one is interested in. Before any sophisticated tests are performed one should check predictions for the energy of the ground state since this quantity is always bounded from below by the exact value of the ground-state energy. Moreover, the energy of the ground state is a quantity, which in systems of few ultra-cold particles can be measured experimentally with high accuracy [31,33].
To test the predictions of the variational method based on this natural quantity we compare the variational energy of the ground state with its counterpart obtained with the exact diagonalization method. The results are presented in Fig. 1, where solid lines represent variational ansatz predictions whereas crosses and squares correspond to the exact-diagonalization predictions (see caption of Fig. 1 for details). Quite obviously, the energy is well reproduced in the limiting cases of g = 0 and g = ∞. However, for the intermediate interactions the energy is clearly overestimated. Moreover, in the perturbation regime of small interactions (g ≈ 0) the slope ∂E(g)/∂g| g=0 is not predicted correctly. These results could suggest that the variational assumption that the ground state of the system can be well approximated with a simple superposition of two many-body eigenstates in limiting cases is maybe too simple. At this point it is worth noticing, that the variational ansatz we use can be essentially improved to make predictions of the ground-state energy much more accurately. The modification is extensively described in [39] and briefly discussed in the Appendix. The improved results obtained in this framework are presented in Fig. 1 by dashed lines. It is clear that the improvement of the resulting energies is essential. Nevertheless, as shown in [39], in this case one loses accuracy in the predictions of the many-body wave functions. Therefore, in further discussion of other quantities the original ansatz Eq. (2) is used and the modified ansatz is adopted only when displaying the energy spectrum in Fig. 1.
Let us note here that also the exact diagonalization method has some problems, mostly in the limit of very strong interactions. It is related to the fact that the resulting energies converge to the exact value very slowly with increasing cutoff N max . Nevertheless, in principle one has a full control on this convergence and can unambiguously indicate a systematic error related to this numerical approximation. However, in cases where convergence is prohibitively slow, the access to a simple ansatz is extremely useful.
Overlap of the many-body ground states.-It is quite natural that in the case of any variational method used to determine the ground state of a many-body problem, a coincidence of energies is not sufficient to claim that the quantum state is predicted correctly. One of the methods to check if the quantum state is reproduced correctly is to calculate its fidelity, i.e., an overlap of the approximate state with the many-body ground state obtained from the exact diagonalization method: Obviously in the case studied, for g = 0 and g → ∞, the fidelity F is equal to 1 since in these limiting cases the wave function is reproduced exactly. For intermediate interactions the fidelity is smaller than 1, and it is presented in Fig. 2. Surprisingly, for equal mass mixture µ = 1 (thin line) the overlap is close to 1 for any interaction, i.e. the wave function of the ground state is reproduced correctly. However, if the mass ratio increases wave function served by the ansatz may have values close to those obtained from the exact method.
To check this hypothesis, in the following we will compare different predictions of the variational approximation with predictions of the exact diagonalization method.
Single particle density.-Apart from the groundstate energy, one of the quantities, which can be measured straightforwardly in experiments, is a spatial density profile of the particles of a given component. Typically it is done by repeating and averaging instantaneous detections of positions of all particles. In principle, in the limit of an infinite number of repetitions, the resulting density approaches the theoretical quantities extracted from the many-body wave function n a (q 1 ) = dq 2 dq 3 dq 4 |Φ g (q 1 , q 2 ; q 3 , q 4 )| 2 , (7a) These profiles can be directly compared with the profiles calculated analogously from the variational ground state of the system |Ψ(g) . Obviously, since the ansatz is based on the proper wave functions in g = 0 and g → ∞, in these limiting cases the predictions of both methods match. If any discrepancies between both predictions exist, one should expect them in the range of interactions where the fidelity F is essentially less than 1. In Fig. 3 we show the density profiles obtained from both methods for g = 2. For equal mass case µ = 1 (left panel in Fig. 3), the exact profile is much flatter than the profile from the variational method. It means that for intermediate interactions variational wave function overestimate contribution from the non-interacting many-body wave Predictions of the variational method (right panels) are consistent with predictions of the exact diagonalization (left panels) for equal-mass µ = 1 (top panels) as well as different mass µ = 10 (bottom panels) systems. However, in the latter case, i.e. for the mass imbalance µ = 10, the probability of finding both fermions in the middle of the trap is overestimated. This fact has consequences in the single-particle profiles (right panel of Fig. 3) where incomplete separation of densities is predicted by variational method. Again, the ansatz results are based on the orginal ansatz as in Fig. 2 and Fig. 3. The positions q1, q2 are measured in natural units of harmonic oscillator /(m b ω), and the two-body density is measured in units of m b ω/ .

function.
When the mass difference between atoms is introduced (right panel in Fig. 3), the density profile of a heavier component is improved. At the same time, the density of a lighter component becomes worse. We checked that this scenario is quite general and it does not depend on statistics, i.e., the result is the same when analogous variational method is adopted for Bose-Bose or Bose-Fermi mixtures.
Although the density profiles predicted by the variational ansatz have some discrepancies when compared to exact results, these differences are rather marginal and should not be of importance in comparison to experimental results. Additionally, we checked that also on the level of a complete single-particle density matrix (not only on its diagonal part) the predictions of the variational ansatz are very close to the exact results. It means that the proposed variational wave function can be safely used to predict any single-particle properties of the system of equal as well as different masses.
Interparticle correlations.-A natural question, which arises at this point, is related to different interparticle correlations that are beyond description of a singleparticle density matrix. Since the ansatz is based on a very simple superposition of two many-body states, in principle, it is not obvious if mutual correlations between particles, which are very sensitive to any change in the many-body wave function, are restored correctly.
To answer this question we concentrate on the simplest two-body correlation, i.e., a two-particle density profile between components defined as Density profiles for interacting system of four fermions of the same and different masses are presented in top and bottom panels of Fig. 4, respectively. As before, the presented results are obtained for the intermediate interactions g = 2, where the fidelity F is essentially lower than 1. It is seen that, in general, the predictions of the variational method are also consistent with exact results. However, some differences are visible, especially for the different-mass systems. Firstly, the variational pair density profiles are much more smeared than the profiles obtained with the exact method. In addition, for the different-mass system, the exact probability of finding both particles in the middle of the trap, in contrast to predictions of the variational method, rapidly drops to zero when the mass ratio µ is increased. This observation is one of discrepancies of the variational ansatz, which may lead to some quantitative differences as compared to experimental data.
Occupations.-One of the less obvious ways of comparing results obtained with different methods is checking the predictions for the occupations of the single-particle orbitals, i.e., the quantities which mathematically are defined for the variational ground state of the system |Ψ(g) as P a (k) = lmn kl; mn|Ψ(g) , For the exact ground state of the system |Φ(g) the definitions are analogous. These quantities are quite interesting in the context of ultra-cold atoms since they can be measured experimentally by an appropriate lowering of the external confinement [28]. Therefore, the theoretical predictions for these quantities can be validated.
In Fig. 5 we present probabilities Eq. (9) calculated for some of the lowest single-particle states as the functions of interactions g for equal (top panel) and different (bottom panel) mass systems. The results based on the variational method (lines) are compared with the probabilities obtained from the exact-diagonalization approach (squares, crosses, etc.). Obviously, in the case of an equal-mass system, both flavors have exactly the FIG. 5. The probabilities Eq. (9) of finding a single fermion in a given single-particle orbital of the harmonic confinement as functions of interactions. In the limit of vanishing interactions the fermions can be found only in the two the lowest orbitals. When the interactions are present other orbitals contribute to the ground state of the system. The predictions of the variational method (solid grey for k = 0, dashed black for k = 1, solid black for k = 2, and dashed grey for k = 3) are roughly consistent with the exact diagonalization results (crosses for k = 0, squares for k = 1, circles for k = 2, and triangles for k = 3). However, for the mass imbalanced systems and stronger interactions the variational method predicts too rapid drop of the ground-orbital contribution below the contribution of the second excited orbital. Here the results are also obtained with the original ansatz as done in Fig. 2, Fig. 3, and Fig. 4. The interaction strength g is measured in units of 3/2 ω 1/2 m −1/2 b . same probabilities. For vanishing interactions particles can be found only in the two lowest states (black solid lines and crosses or black dashed lines and squares for states with k = 0 or k = 1, respectively). As interaction increases, both probabilities decrease and higher singleparticle states become partially occupied. In this case, predictions of the variational method, although not perfect, reproduce results from the exact method quite well.
The situation changes significantly, when a mixture of different masses is considered. In this case, predictions of both methods are roughly consistent only for the heavier component. For the lighter component, the occupation of the lowest single-particle orbital rapidly drops with the increase of the interactions and the third orbital becomes significantly occupied. Moreover, at some moment the occupation of the ground orbital becomes less probable than the occupation of the third state. This behavior of probabilities for lighter component is predicted by the variational ansatz. However, for small interactions a mentioned drop is too slow, whereas for stronger interactions (around g = 2) it is too rapid. Nevertheless, differences between exact diagonalization predictions and variational approach are not essential. It means that also in the case of different-mass systems the variational ansatz can still be used for qualitative predictions when occupations of different single-particle orbitals are considered.
This extremely good agreement between the predictions of the variational ansatz and the exactdiagonalization approach in the case of the heavy component is related to the fact that in both limiting cases (g = 0 and g → ∞) the heavy particles are located in the middle of the trap. The situation is different for light particles, i.e., in the limit of strong repulsion the light particles are pushed out from the middle of the trap. This implies that for the light particles interactions have a stronger effect spatially. As a consequence, it is considerably more difficult to capture this effect by the ansatz constructed as a superposition of the limiting wave functions. This effect is directly reflected in the occupations of the single-particle orbitals.

IV. FINAL REMARKS
Other statistics.-Although the results presented are related to the fermionic mixtures, to obtain a wider perspective on the problem of accuracy of the variational ansatz, it is worth considering different kinds of mixtures of four quantum particles. Both methods, i.e., the variational ansatz and the exact diagonalization approach, can be easily adopted for mixtures of two kinds of bosons or one kind of boson and one kind of fermion. Formally, the only difference has to be introduced in (anti)symmetrization definitions in Eq. (5). Of course, these changes may have (and typically do have) decisive consequences for the results obtained.
We have performed appropriate calculations for Bose-Bose and Bose-Fermi mixtures under the assumption that the bosons within a given flavor do not interact and the only non-vanishing interaction is present between different components. This assumption gives us a simple comprehensive tool for testing the role of quantum statistics. The strongly interacting states for the equal mass case in Bose-Fermi mixtures have been a subject of several recent discussions [54][55][56]. We note that in Bose-Bose mixtures with no interactions within a given flavor, the strongly interacting wave function cannot be found by building it on a basis of totally antisymmetric wave function, but other techniques to obtain it have been discussed recently [43][44][45].
While we do not present the full results of our calculations here, some of the results obtained in this way were already mentioned previously. From our numerical tests and comparisons some general conclusions about the role of the statistics can be given. Independently of the statistics, the simple variational ansatz works surprisingly well and it can be safely used for simple qualitative and quantitative predictions when single-particle observables are considered. In fact, the case that we have presented here with two kinds of fermions is that in which the comparison between the numerically exact and the variational method is the worst. For other compositions of the particles the variational method agrees even better with the numerical results. We do caution though that whenever higher interparticle correlations are considered, one should be very careful since the predictions of the variational method proposed can be overestimated.
Improving the ansatz.-It is quite obvious that in principle the variational probe function (2) could be extended by superposing additional many-body statefor example the many-body ground state obtained numerically for the interaction g for which the accuracy is the worst. Although such extension is possible it requires some numerical effort to obtain an additional many-body state. Therefore, a lot of the beauty and simplicity of the idea may be quickly lost. Nevertheless, this direction would be necessary if larger number of particles were considered.
Conclusions.-In this paper we compared predictions of the interpolatory ansatz introduced in [39] with the numerically exact method of diagonalization of the many-body Hamiltonian. Surprisingly, the simple assumption that the ground state of four interacting fermions of different masses can be well approximated by a superposition of two many-body ground states obtained in the limits of very strong and vanishing interactions, is sufficient to describe many different properties of the system. Obviously, in this simplified description some discrepancies are present for the intermediate interactions, but they are rather small and not decisive in the view of a quite drastic simplification.
As a consequence, the optimized energy reduces to The above estimation of the energy turns out to be insufficient in the limit of strong as well as weak interactions [39]. It can be improved by requiring the correct slope of the energy as a function of −1/g. It is shown that up to the first-order expansion the slope of the energy in the strong interactions regime can be calculated exactly and it is given as [7]: where K 0 = Ψ 0 |H ab |Ψ 0 /g. For a given K ∞ opt we can now find a new value of Ψ 0 |Ψ ∞ , which can be inserted in the expression for the optimized energy (11). In this way one obtains much better estimation of the groundstate energy. Even though the modified ansatz reproduces the energy much better than the original ansatz (see Fig. 1), it comes with the cost that the ground-state wave function is no longer known. As a consequence, the modified ansatz is useful only for performing a better estimation of the energy.