Pair-correlation ansatz for the ground state of interacting bosons in an arbitrary one-dimensional potential

We derive and describe a very accurate variational scheme for the ground state of the system of a few ultra-cold bosons confined in one-dimensional traps of arbitrary shapes. It is based on assumption that all inter-particle correlations have two-body nature. By construction, the proposed ansatz is exact in the noninteracting limit, exactly encodes boundary conditions forced by contact interactions, and gives full control on accuracy in the limit of infinite repulsions. We show its efficiency in a whole range of intermediate interactions for different external potentials. Our results manifest that for generic non-parabolic potentials mutual correlations forced by interactions cannot be captured by distance-dependent functions.

At the moment the quantum mechanics was rigorously formulated, it became obvious that understanding of different properties of quantum many-body systems is incomparably more challenging than their one-and twoparticle counterparts 1 . The main source of this obstacle is of course the tensor product structure of the manybody Hilbert space growing exponentially with the number of particles. That simple fact, already for systems containing only several particles in one spatial dimension, leads directly to an impassable fiasco of any accurate numerical approach willing to take into account all possible many-body configurations coupled by mutual interactions. In fact, even after firm reduction of allowed states (justified only for relatively weak interactions), direct methods like the exact diagonalization 2 are generally time-consuming and not efficient. Of course, some clever modifications of these approaches may shift the boundary of their applicability [3][4][5][6][7] . Still, however, the main stumbling block, i.e., the exponential growth of the Hilbert space, remains unchanged. Moreover, only few exactly solvable models for many-body systems are known like: the Moshinsky model 8,9 , the Lieb-Liniger model 10,11 , the McGuire-Gaudin-Yang model [12][13][14][15][16] , the Calogero-Sutherland model 17,18 , the Yang-Baxter integrable models 19 , or the one-dimensional model with some particular combination of contact and long-range interactions found recently 20 . In this context, broad research on effective numerical methods and techniques dedicated to different many-body systems are continuously ongoing [21][22][23][24][25] .
Besides all sophisticated methods dedicated to quantum many-body systems, there is one strategy that is well-suited to the problem of finding isolated many-body ground states. The method is very natural since it is based on intuitive variational arguments 26 . Simply, one needs to postulate a family of trail functions depending on a certain set of variational parameters. These parameters are optimized (analytically or numerically) to grant possibly minimal expectation value of the many-body Hamiltonian. If the trial family is luckily chosen appropriately, then the obtained wave function approximates the actual ground-state of the system very closely. Thus, all the ground-state properties are easily accessible with relatively low computation costs. Up to now, different approaches in this fashion have played a huge role in investigating different properties of various quantum systems.
It is worth emphasizing that direct insight into the mathematical structure of inter-particle correlations served by variational approaches should be recognized as a huge advantage having crucial importance. This is true especially if the state-of-the-art experiments with several ultra-cold atomic systems are considered 27,28 . In this area, different experimental techniques enable one not only to precisely control and tune different parameters of a system but also to measure high-order correlations between particles with tremendous accuracy 29  www.nature.com/scientificreports/ Therefore, having well-justified theoretical predictions for the state's structure is a mandatory requirement for a clear understanding of the system properties [30][31][32] .
In the context of variational approaches applied to few-body problems, some effort was devoted recently to construct an appropriate variational trial wave function for particles confined in traps being close to parabolic shape [33][34][35][36][37][38] . Starting from the well-known exact solution for two interacting particles 39 , following the Jastrow idea 40 , there were proposed reasonable many-body wave functions for bosons 41,42 as well as fermions 43 . Some sort of generalization to multi-component systems was also given 42,44 . All these attempts have one fundamental obstacle -they cannot be straightforwardly adapted to other confinements since then they inappropriately capture the behavior of the system for strong interactions. In this work, by a deep investigation of a few-boson system confined in the external potential of an arbitrary shape, we show how to overcome this restriction. Notably, first, we elaborate on a two-particle problem and we define the ansatz that is exact in the limit of vanishing and infinite repulsions, while for intermediate interactions it nicely reproduces the actual ground-state of the system. Based on this experience, we show the strategy for obtaining the accurate ansatz for a larger number of particles. We demonstrate its effectiveness on different model systems with single-and double-well potentials over the entire interaction regime.
At this point, we want also to point out that the method presented here for confinements being far from the parabolic is not the only possible path to tackle the problem. For example, the problem of impurities confined in an asymmetric double-well potential was considered recently in the framework of the variational interpolatory ansatz 45 .

Formulation of the problem
We consider a closed system of N indistinguishable spinless bosons of mass m, interacting mutually via δ-like contact interactions, and confined in a quasi-one-dimensional external trap. The Hamiltonian of the system has a form: where V(x) determines a shape of the external potential and g controls the strength of interactions. Experimentally, such systems are attainable with ultra-cold dilute atomic gases 32 . A quasi-one-dimensional confinement is obtained by applying very strong confinement in two perpendicular directions and, besides utilization of the Feshbach resonances, may be used to control effective strength of interactions 46 .
In general, we aim to find the ground-state wave function of the system for possible general confinement V(x). Although our strategy can be applied almost to any reasonable confinement, motivated by recent experimental implementations and theoretical proposals [47][48][49][50][51][52][53][54] , in this work we focus on a quite general class of double-well potentials. They are modeled by a combination of a harmonic potential, a gaussian barrier in the middle, and a gradient field having a linear slope. It can be written as: By controlling three independent parameters , δ , and ξ one can tune the trap to the pure harmonic oscillator ( = 0 ), symmetric ( ξ = 0 ) and non-symmetric ( ξ = 0 ) double-well potential of controlled barrier height and width δ > 0 . To make further analysis clear, in the following we express all quantities in terms of natural units of the pure harmonic oscillator, i.e., energy, length, and interaction strength are expressed in units of ℏω , √ ℏ/mω , and ℏ 3 ω/m , respectively. In this convention the parameters , δ , and ξ are dimensionless (they are expressed in units of energy, length, and force, respectively). For convenience, but without losing generality, we set the energy shift V 0 to the value ensuring the condition V (0) = 0.
Before we present our variational scheme let us note here that even if the corresponding single-particle Hamiltonian is diagonalized and its eigenorbitals ϕ i (x) are known, the exact solution of this many-body problem is known only in a few specific cases. First, in the absence of interactions ( g = 0 ) all bosons occupy the lowest single-particle orbital ϕ 0 (x) and thus the ground-state wave function reads �(x 1 , . . . , x N ) = i ϕ 0 (x i ) . Second, for two bosons confined in a harmonic trap and arbitrary interactions the known solution has a form �(x 1 , is the confluent hypergeometric function. The parameter ǫ g is determined (via transcendental relation) by the interaction strength g. The detailed derivation of this solution can be found in the original paper by Busch et al. 55 or in a more pedagogical work 56 . Finally, as shown by Girardeau 57 , in the limit of infinitely strong repulsions ( g → +∞ ), the N-boson ground-state wave function is also known exactly. Indeed, it can be directly obtained from the wave function of N non-interacting fermions confined in the same external potential by applying an appropriate symmetrization. In all other cases, one can find the ground-state wave function only in an approximate fashion. For small or very strong interaction strengths it can be done perturbatively by expansion around known exact solutions or by performing numerical calculations based on the exact diagonalization. The most challenging (and therefore the most interesting) are of course cases with intermediate interactions where any expansion to any order becomes insufficient. In consequence, any numerical method is not reasonably well converged.

Variational ansatz
Here we propose a relatively simple but in many cases very accurate variational ansatz for the bosonic wave function. It is based on the conviction that the dominant part of inter-particle correlations can be expressed by two-body correlations only. Since the bosonic wave function must be symmetric under the exchange of any two particles, this assumption leads directly to the most general form of the wave function: where ϕ(x) and �(x, y) are some unknown single-and two-particle orbitals describing the state of the system and N is a numerical constant assuring appropriate normalization. This form of the many-body wave function can be viewed as a specific generalization of the celebrated Jastrow ansatz 40 adapted to confined systems. Its specific form was used previously in the case of harmonic confinement for bosons as well as for fermions 42,43 . It can be shown straightforwardly 10 that the contact interaction part of the Hamiltonian (1) maybe directly included to the form of the wave function. Indeed, if the pair correlation function �(x, y) fulfills the condition then it is sufficient that the wave function � G (x 1 , . . . , x N ) is found for the ground state of the non-interacting system, i.e., the system described by the Hamiltonian Strictly speaking, if we build a variational family fulfilling the condition (4) then the variational approximation of the ground-state wave function is found by minimizing the non-interacting energy functional To propose the variational family as accurate as possible, first we rewrite the two-particle Jastrow orbital �(x, y) to a more convenient form. We express it in terms of yet unknown antisymmetric function D(x, y) and a single real parameter α as where By introducing above redefinition we automatically assure that the pair-correlation function �(x, y) fulfils contact condition (4) for arbitrary value of α and for arbitrary function D(x, y) which is antisymmetric under exchange of variables, D(x, y) = −D(y, x) . In this way we included automatically all conditions forced by mutual interactions. Thus, we can use the family (after probing different functions D(x, y)) to minimise the energy functional (6).
Although one has very wide freedom in choosing the antisymmetric function D(x, y) and the single-particle orbital ϕ(x) , it is very convenient to tailor them in such a way that known exact results for particular systems are well-captured by the ansatz. Therefore, now we will impose conditions on the trial family forced by limiting cases.
The simplest condition is imposed by considering the non-interacting limit ( g → 0 ). We know that in this limit all bosons occupy the lowest orbital ϕ 0 (x) of a corresponding single-particle problem. Thus, the variational wave function (3) captures this limit appropriately provided that the single-particle orbital ϕ(x) is chosen as the lowest orbital ϕ 0 (x) . Indeed, in the limit g → 0 the variational parameter α → +∞ and therefore the two-particle orbital becomes trivial, �(x, y) = 1.
Relatively more challenging is to fulfill conditions originating in the limit of infinite repulsions ( g → +∞ ). It is clear that this limit is quite well mimicked by the ansatz for vanishing α . Indeed, for α → 0 one finds � α (x, y) → 1 and therefore (due to the expansion e −αz ≈ 1 − αz) We know however that in this limit the exact form of the ground-state wave function can be obtained by exact mapping from the ground state of non-interacting fermions. It means that the exact bosonic many-body wave function is expressed as a product of the Slater determinant formed with the lowest single-particle orbitals and a chain of sign functions 57 : www.nature.com/scientificreports/ Consequently, the pair-correlation function D(x, y) should be chosen such that the function (8) may reconstruct the exact solution (9) as close as possible. Of course, there is no general prescription how to construct function D(x, y) for arbitrary confinement. However, in the case of general double-well traps (2) studied here some rigorous hints can be formulated. The first hint clarifies when a pure harmonic trap potential is considered. Then the exact wave function has a simple form Therefore, for harmonic confinement, by taking correlation function in a very simple form of D(x, y) = x − y one reconstructs with (8) the exact solution rigorously. This kind of pair-correlation function was exploited in previous works 42,43 . The second hint can be found by considering the complementary system of two particles in the confinement of arbitrary shape. Then, the exact solution in the limit of infinite repulsions is expressed by two the lowest single-particle orbitals ϕ 0 (x) and ϕ 1 (x) as: where π(x) = ϕ 1 (x)/ϕ 0 (x) . It means that by choosing D(x, y) = π(x) − π(y) in our variational scheme one reconstructs solution of the two-body problem in arbitrary trap for infinite repulsions rigorously. At this step it is also worth to consider further improvement for finite interactions which does not change applicability in the infinite repulsions limit. Namely, one can extend proposed correlation function by additional variational parameter β responsible for an effective rescaling of single-particle orbitals and define In this way, the ansatz is more flexible and may give improved results for intermediate interactions.
It is worth noticing that in the case of pure harmonic confinement π(x) = x . It means that the parameter β is redundant with rescaling α and therefore it becomes irrelevant. This all means that the proposed form of the correlation function (12) defines very reasonable yet not very complicated ansatz for N = 2 particles in almost any trap. In the following section, we discuss its advantages over the naive ansatz proposed previously based on a harmoniclike substitution π(x) = x.
In the case of a larger number of particles, one needs to take one step further in defining the variational ansatz and to release constraints for function π(x) since it is no longer related to the ratio of the lowest two eigenfunctions of the single-particle problem. This step is essential since for any non-harmonic confinement and a large enough number of particles, the function π(x) must gain substantial corrections from higher single-particle orbitals occupied by particles. Now we will show how to generalize the correlation function D β (x, y) to make it more flexible for different confinements but keeping its appropriateness in limiting cases discussed above.
Since in our work we consider only smooth and bounded potentials it is clear that function π(x) can be expanded in Taylor series π(x) = k p k x k . It immediately means that, up to irrelevant constant factor p 1 β , the correlation function has a form where P k (x, y) is (k − 1)-th order polynomial of the property (x − y)P k (x, y) = x k − y k . Note that for any mirror-symmetric external potential ( ξ = 0 ) the function π(x) is odd ( π(−x) = π(x) ). Then, only odd-k terms are present in its decomposition and therefore the correlation function D β (x, y) is also simplified.
From the obtained structure of the correlation function originating in a two-orbital approach (13) we deduce possible paths of generalization. First, more accurate but numerically very challenging and numerical-timeconsuming, relay on converting decomposition coefficients p k to independent variational parameters. In this approach, the final ansatz of n-th order is defined as where b k is k-th element of variational parameters vector b.
The second possibility, which we widely exploit in the following, goes in opposite direction. Namely, we associate all expansion coefficient b k with a single variational parameter β (we set b k = β k ) but we release a rigid polynomial structure of the expansion (14). Consequently we write the correlation function as www.nature.com/scientificreports/ where the variational function F β (x, y) obeys four conditions: (i) it is symmetric under exchange of variables, F β (x, y) = F β (y, x) ; (ii) it saturates on unity for vanishing variational parameter β , F 0 (x, y) = 1 ; (iii) when Taylorexpanded in β , consecutive terms should be expressed as polynomials of x and y of consecutive powers; (iv) for mirror-symmetric potentials ( ξ = 0 ) only even polynomials contribute to the function F β (x, y) , i.e., terms with even powers of β only. The conditions (iii-iv) are direct consequences of the explicit form of the expansion (13). These conditions mean that proper construction of the correlation function F β (x, y) should treat even and odd polynomial expressions as independent rather than related -intensity of the odd part is controlled mostly by an asymmetry of external potential. Now we are ready to propose the variational family of trial functions obeying desired conditions. In our work, we postulate one of the simplest forms which can be applied for any external potential having the form (2). It is constructed from two the lowest even polynomials and supplemented by the lowest odd polynomial to capture breaking of the mirror symmetry for cases with ξ = 0 . The correlation function has a form In this definition, mentioned independence of even and odd polynomials is provided by the additional parameter γ which vanishes for mirror-symmetric scenarios. In general, it can be treated as an additional variational parameter. However, we checked that for a given number of particles and a given potential shape, it is sufficient to establish its value in the limiting case of infinite repulsions and then use it as the fixed parameter for all the other interactions. Of course, in the case of N = 2 particles, it is not needed to use the ansatz (16) since the previous approach based on the definition D β (x, y) = π(βx) − π(βy) gives similarly good results and is significantly much simpler in numerical adaptation. For more complicated confinements, like triple-or more-well potentials, one should try to extend the function D β (x, y) by higher-rank polynomial expressions in the exponent, keeping all the conditions for function F β (x, y) always fulfilled.

Validation method
To show that the proposed variational ansatz (16) describes adequately properties of the many-body system confined in different traps we compare its predictions to the results obtained numerically by the exact diagonalization of the many-body Hamiltonian. For this purpose we use the optimised exact diagonalization method with the harmonic oscillator basis 6 . Since we deal with external confinements having parabolic shape far from the center, this approach assures that we get states having appropriate asymptotic behaviour for x → ±∞ . For a given number of particles N, interaction strength g, and parameters of the external confinement, the many-body Fock basis is constructed as follows. First we form a set of K + 1 single-particle orbitals {ϕ where H k (x) is the Hermite polynomial. Note that plays a role of an optimalization parameter and is not dependent on frequency ω in (2). Then we construct the Fock basis over these states. Namely, in the secondquantized occupation notation the basis is formed by Fock states |F ℓ � = |n 0 , n 1 , . . . , n K , 0, . . .� (the index ℓ enumerates Fock states) with the cut-off condition K k=0 k n k < K . In the next step, all matrix elements of the Hamiltonian H ℓℓ ′ = �F ℓ |H|F ℓ ′ � are calculated and the corresponding matrix is diagonalized. In this way approximate ground-state energy E { } 0 and associated eigenvector representing the many-body ground state |G { } � are found. Finally, we repeat all these calculations for different values of the optimalization parameter to find the lowest ground-state energy E 0 for a given cut-off parameter K. By increasing K, the final results are systematically improved and they converge to their ultimate (physical) values.
In our comparison we focus on systems with N ∈ {2, 3, 4} particles and the most natural quantity -the singleparticle reduced density matrix defined in the position representtion as This quantity encodes all possible results which can be obtained in arbitrary single-particle measurements. The diagonal part corresponds to the density profile of atomic cloud, while the off-diagonal elements of the matrix are responsible for spatial coherence of the system on a single-particle level. They can be quantified in many different ways. To deliver a reasonably complete comparison, in our work we consider three different cuts of the reduced density matrix: the diagonal, the anti-diagonal, and the horizontal. They are defined as By comparing these three distributions obtained via the ansatz and the exact numerical calculations we can give a quite comprehensive validation test.

Two-particle case
We start validation of the ansatz in the simplest case of N = 2 particles. As was mentioned previously, in this case it is not necessary to use the general ansatz (16) but it is sufficient to utilize its simpler form D β (x, y) = π(βx) − π(βy) (see Eq. (12)). In Fig. 1, we plot different cuts of the single-particle density matrix ( ρ(x) , m(x), and h(x) on consecutive columns) predicted by the ansatz (different symbols) and we compare them www.nature.com/scientificreports/ with corresponding results obtained via the exact diagonalization (solid lines). To make a comparison comprehensive we analyze three different interaction strengths and two different double-well-like confinements (first and second row, respectively). It is clear that the ansatz correctly predicts shapes of all three distributions in a whole range of interactions not only in the case of symmetric double-well confinement ( ξ = 0 ) but also in the case of a very asymmetric scenario ( ξ = 0.2 ). For comparison, in the last column in Fig. 1 we display predictions for the density distribution ρ(x) of the same ansatz but with naive substitution D β (x, y) = x − y which was shown to be perfectly suited in the case of confinements close to the harmonic one 42,43 . It is rather obvious that the naive ansatz, even in the case of symmetric double-well confinement, is not optimal. Especially this is the case for positions being close to the double-well minima. The situation is even worse for asymmetric traps where this simplified ansatz is completely unreliable. This evident inappropriateness for asymmetric cases comes mainly from the fact that here the correlation function D β (x, y) shouldn't be mirror-symmetric. Mentioned breaking of mirror-symmetry is directly encoded in function π(x) which is built from two of the lowest single-particle orbitals, while it is completely neglected by setting D β (x, y) = x − y.
Appropriateness of the ansatz studied and clear deviations of the naive approach can be visualized further when, instead of different cuts of the reduced single-particle density matrix, its global properties describing inter-particle correlations are considered. They are nicely quantified by the entanglement entropy S which is one of natural measures of non-classical correlations between particles. Particularly, in the case of two indistinguishable particles, the general criterion of their entanglement based on the entanglement entropy was formulated 58 . It can be easily calculated if the decomposition of the reduced density matrix to its natural orbitals {u i (x)} is known. Namely, in the basis of natural orbitals, the density matrix ρ(x, x ′ ) is expressed as a sum of projectors where each r i is the eigenvalue of ρ(x, x ′ ) associated with the eigenorbital u i (x) . Then, having eigenvalues {r i } in hand, one defines the von Neumann entanglement entropy S = − i r i ln r i . For the two-particle scenario considered here, we compute the eigenvalues {r i } directly from the wave function, i.e., we diagonalize the matrix � ij = [�x� G (r i , r j )] on a dense grid x = r i+1 − r i and we find a set of its eigenvalues {k i } . Then the eigenvalues of ρ(x, x ′ ) are simply given as i = k 2 i (see for example 59 ). In Fig. 2 we display the entanglement entropy S as a function of mutual interactions g for three different traps (displayed in the inset). We compare results obtained with two different approaches: the exact diagonalization (solid black lines) and the proposed variational ansatz (red symbols). Both approaches give almost the same results for a whole range of interactions and any external confinement considered. For comparison, in the case of shallow and asymmetric double-well confinement (trap B) we also display results obtained with the naive approach based on substitution D β (x, y) = x − y (black crosses). In this case, deviations from the exact results www.nature.com/scientificreports/ and substantial underestimation of inter-particle correlations are clearly visible. This result supports previous findings based on the single-particle density profile.

Larger number of particles
As argued before, in the case of a larger number of particles one needs to release stiffness of the ansatz and use a slightly more general family defined according to (16). To show that the ansatz is very accurate also in these cases we perform a detailed analysis of systems containing N = 3 and N = 4 particles confined in three very different confinements being far from the parabolic trap: symmetric double-well potential, asymmetric double-well potential, and essentially asymmetric single-well trap (detailed parametrization is given in the caption of Fig. 3). In all these cases we investigate the accuracy of the ansatz for interactions for which still we are able to perform exact diagonalization of the Hamiltonian and thus establish an appropriate benchmark. Additionally, we also check predictions of the ansatz in the limit of infinite repulsions ( g → ∞ ) where the exact form of the many-body ground state is known. First, in Fig. 3, we display predictions for the diagonal part of the single-particle density matrix, i.e., the density profile n(x). It is clear that the ansatz's predictions are very close to the exact results for all three types of traps and any interaction strength considered. Note that along with increasing interactions the density distribution undergoes substantial transformations. This is especially the case for asymmetric traps where the density transits from a very asymmetric distribution (for weak interactions) to much more regular (for very strong interactions). Importantly, this behavior is perfectly captured by the ansatz. This very nice property of the variational scheme proposed is a direct consequence of its construction which by definition is tailored to reconstruct appropriately the limit of infinite repulsions. At this point let us also mention that obtaining numerically exact results via direct diagonalization of the Hamiltonian becomes numerically challenging when the number of particles is increased. Therefore, in the case N = 4 , we compare results for relatively weaker strengths than in the case of N = 3 particles. Contrary, the numerical complexity of the ansatz approach is almost insensitive to the interaction strength and the number of particles. Thus, it can be successfully used for cases being far beyond the feasibility of the numerical diagonalization approach. Quality of the ansatz's predictions is also granted when instead of the diagonal part, different off-diagonal properties of the single-particle reduced density matrix are considered. We checked that two other quantities we focus on in our report, namely m(x) and h(x), are perfectly reconstructed by the ansatz for any shape of the external trap and any interaction strength. For example, in Fig. 4, we show corresponding results for the most challenging confinement -asymmetric double-well trap (middle case in Fig. 3). For this trapping configuration, the exact diagonalization approach is on the border of feasibility and obtaining well-converged results requires huge numerical and computational resources. In contrast, the numerical complexity of the approach based on the ansatz (16) does not substantially depend on the external trapping and even in this demanding case is easily manageable.
At this point it is worth admitting that in contrast to other quantities, obtaining the von Neumann entanglement entropy S from the ansatz is not an easy task. By the construction, the ansatz serves the ground-state as a wave function in the position domain. Thus, it is very easy to obtain all reduced quantities also in this  www.nature.com/scientificreports/   www.nature.com/scientificreports/ representation. Entanglement entropy is not of this kind since it is obtained after diagonalization of the density matrix which in the position domain requires numerical integrations on a very dense grid. Thus, for systems with more than N = 2 particles obtaining this particular quantity from the ansatz is not efficient.

The ansatz as a convenient validator
Up to now, our validation method was based on comparison with the results obtained via numerically exact diagonalization of the Hamiltonian. Of course, this method has limited capabilities since the diagonalization is always performed on a cropped basis formed by states with the lowest excitations. Thus it is not credible when highly-excited states have significant participation to the ground state. This is one of the fundamental reasons why in our previous comparisons we considered only weak enough interaction strengths for which the method is wellconverged, i.e., further inclusion of highly excited states does not change obtained results. To expose this issue, in Fig. 5, we compare the ground-state energy as predicted by the ansatz and the exact diagonalization method for different interactions and different numbers of particles. It is clear that along with an increasing number of single-particle orbitals K from which the Fock basis is built, successive ground-state energies predicted by the exact diagonalization converge to their true values. At the same time, horizontal lines represent corresponding ground-state energies predicted by the ansatz. It is evident that typically the ansatz approximates the ground state much better than the diagonalization approach. Only for a small number of particles, relatively small interactions, and large enough Fock basis the ansatz estimates the ground-state energy slightly worse than the numerical diagonalization. It means that the ansatz proposed is not only much more convenient and faster than a direct diagonalization but in most cases gives results being closer to the physical reality. From this point of view, the ansatz should be recognized as an appropriate tool for the validation of other methods rather than the opposite.

Conclusions
In our work, we have shown how to get a very precise approximation of the many-body ground state of a onedimensional system of several ultra-cold bosons in terms of the pair-correlated variational scheme. By detailed analysis of the contact interaction terms and a role played by external confinement, we give a straightforward recipe for appropriate construction of the variational family. Additionally, we have pointed out sources of potential discrepancies and we explained how to control them by including different variational terms. In principle, for a given external potential, the main idea is to treat the exact solution for infinite repulsions as a conductor benchmarking the accuracy of a whole method. The effectiveness of the approach proposed enabled us to conclude that the shape of the trapping potential strongly determines the spatial structure of particles' correlations, including the degree of entanglement. As consequence, we argued that other, recently used approaches 39,42,43 , assuming pair-correlators as simple distance-dependent functions are insufficient for non-harmonic potentials. From this point of view, the presented results broaden the recent discussion on mutual correlations in systems of few ultra-cold atoms and open paths for additional searching for effective and accurate approximate methods. Although our work is devoted to the simplest problem of spinless bosons, the method presented can be generalized to the case of a larger number of components and also having different statistics, similarly as it was done previously for the naive ansatz with D β (x, y) = x − y 42 . It can be done as long as the system considered is one-dimensional and mutual interactions have zero-range nature. Note that along with increasing cut-off K the numerically exact results converge from above to the ansatz's predictions. This indicates that the variational scheme proposed, although relatively simple and numerically less demanding, typically gives very accurate and better results than very heavy calculations in the framework of the matrix diagonalization. In the presented examples, calculations are performed for the same trap as in Fig. 4 (corresponding to the second column in Fig. 3).