Sorting Fermionization from Crystallization in Many-Boson Wavefunctions

Fermionization is what happens to the state of strongly interacting repulsive bosons interacting with contact interactions in one spatial dimension. Crystallization is what happens for sufficiently strongly interacting repulsive bosons with dipolar interactions in one spatial dimension. Crystallization and fermionization resemble each other: in both cases – due to their repulsion – the bosons try to minimize their spatial overlap. We trace these two hallmark phases of strongly correlated one-dimensional bosonic systems by exploring their ground state properties using the one- and two-body density matrix. We solve the N-body Schrödinger equation accurately and from first principles using the multiconfigurational time-dependent Hartree for bosons (MCTDHB) and for fermions (MCTDHF) methods. Using the one- and two-body density, fermionization can be distinguished from crystallization in position space. For N interacting bosons, a splitting into an N-fold pattern in the one-body and two-body density is a unique feature of both, fermionization and crystallization. We demonstrate that this splitting is incomplete for fermionized bosons and restricted by the confinement potential. This incomplete splitting is a consequence of the convergence of the energy in the limit of infinite repulsion and is in agreement with complementary results that we obtain for fermions using MCTDHF. For crystalline bosons, in contrast, the splitting is complete: the interaction energy is capable of overcoming the confinement potential. Our results suggest that the spreading of the density as a function of the dipolar interaction strength diverges as a power law. We describe how to distinguish fermionization from crystallization experimentally from measurements of the one- and two-body density.


Hamiltonian, One-and Two-Body Density
In order to discuss the stationary properties of the ground state (GS) of crystalline and fermionized bosons, we consider the time-independent many-body Schrödinger equation, Here, |Ψ〉 is the many-body ground state, E its energy, and Ĥ the N-particle Hamiltonian in dimensionless units 56 , to be the external harmonic trap. The term W(x i − x j ) is the interaction potential. All quantities are dimensionless and expressed in harmonic oscillator units. To ensure that the system is in the quasi-1D regime, we assume strong confinement in the transversal direction, providing a cigar-shaped atomic density. The contact interactions read, where λ is the interaction strength determined by the scattering length a s and the transverse confinement frequencies 57 . For long-ranged dipolar interactions, we have where g d is the dipolar interaction strength and α is a short-range cut-off to avoid the divergence at x i = x j . Repulsive interactions could be obtained in a quasi-1D BEC by imposing an external magnetic field to align all the dipole moments of the atoms 18 . This simple approximation to the one-dimensional dipole-dipole interaction potential in Eq. (4) is justified for the moderate to large interaction strengths and large inter-particle distance with respect to the harmonic-length of the transversal confinement 14,15,18,19,58 , that we focus on in the present work. For such interaction strengths the dipole-dipole interaction potential is well-approximated by the |x i − x j | −3 tail in Eq. (4), see 59 . Moreover, we have verified the consistency of the approximation in Eq. (4) for the same choice of cutoff parameter, α = 0.05, by a direct comparison to a dipole-dipole interaction augmented with an additional contact interaction potential, see Ref. 34 . A rigorous discussion of the dipole-dipole interaction potential in one and two spatial dimensions can be found in Ref. 59 . Here, for the sake of simplicity, we will focus on quasi-one-dimensional systems, with N = 4 interacting bosons for all our calculations and consider repulsive interactions, i.e., λ > 0 and g d > 0, exclusively. In the following we discuss the reduced one-body density matrix, defined aŝ † (1) Its diagonal, is simply the one-body density. As a precursor of correlation effects that may be present in the state |Ψ〉 of the system, we use the eigenvalues ρ i (NO) of the reduced one-body density matrix ρ (1) in Eq. (5). For this purpose, we write ρ (1) in its eigenbasis: The eigenvalues ρ i (NO) and eigenfunctions Φ i (x) are referred to as natural occupations and natural orbitals, respectively. If only a single eigenvalue ρ i (NO) is macroscopic, then the state |Ψ〉 describes a Bose-Einstein condensate 60 . The case when multiple eigenvalues ρ i (NO) are comparable to the number of particles N is referred to as fragmentation 43,44,53,54,61-63 . In the following, we will also use the two-body density ρ (2) to characterize crystallization and fermionization. It is defined as (2) 1 2 1 2 1 2 The two-body density quantifies the probability to detect two particles at positions x 1 and x 2 .

numerical Method
The computation of the exact many-body wave function is a difficult problem. To attack the many-body Schrödinger equation, Eq. (1), we use the time-dependent Schrödinger equation, t with a Wick's rotation t → −iτ, i.e., a propagation with imaginary time. We expand the many-body wavefunction |Ψ〉 of N interacting bosons in a complete set of time-dependent permanents | → 〉 = | ... 〉 n t n n t ; , , ; M 1 with at most M single-particle states or orbitals. The MCTDHB ansatz for the many-body wave function is thus Here, the permanents | → 〉 n t ; are symmetrized bosonic many-body states that are also referred to as "configurations". The sum in Eq. (10) runs on all configurations → n of N particles in M orbitals. The number of permanents and In second quantized representation the permanents are given aŝ is the bosonic creation operator which creates a boson in the time-dependent single particle state φ x t ( , ) k . Equation (10) spans the full N-body Hilbert space in the limit of M → ∞. For practical computations, we restrict the number of orbitals and require the convergence of our observables, like the one-and two-body density matrix, with respect to the number of single-particle states M.
A set of coupled equations of motion for, both, the time-dependent expansion coefficients → C t ( ) and φ x t ( , ) k . Using MCTDHB, both, the coefficients and orbitals are variationally optimized 64 . MCTDHB is thus fundamentally different from exact diagonalization, i.e., an ansatz built with time-independent orbitals. It can be demonstrated that MCTDHB delivers solutions of the Schrödinger equation at a significantly increased accuracy in comparison to exact diagonalization approaches when the same number of single-particle basis states is employed, see Refs. 53,65 . for a demonstration with the harmonic interaction model and Appendix A for a demonstration with dipole-dipole interactions, i.e., the Hamiltonian in Eqs. (2) and (4). Despite the accuracy of MCTDHB for weakly interacting particles, for strong interactions a large number of orbitals is required to describe the system accurately. In the case λ → ∞, the Bose-Fermi mapping provides analytical solutions that can be compared to the numerical results. These difficulties to converge MCTDHB results for strong contact interactions may be an instance of the discussion provided in Ref. 66 .
We solve the set of coupled MCTDHB equations using the MCTDH-X software 46,[52][53][54] . When imaginary time is used, the propagation of an initial guess function converges to the ground state of the system, and the stationary properties of the system can be investigated.

fermionization vs crystallization
We now discuss our findings on the fermionized and the crystalline state of parabolically trapped one-dimensional ultracold bosons. We first independently characterize fermionization and crystallization from a "many-body point of view", see Sec. 4.1 and Sec. 4.2, respectively. Thereafter, we investigate how to sort the one, fermionization, from the other, crystallization, in Sec. 4.3. Here and in the following, we used the term "many-body point of view" to highlight that our considerations go beyond an effective single-particle or mean-field description of the state. fermionization. Bosons fermionize when they feel an infinitely repulsive contact interaction in one spatial dimension. For fermionized bosons, the total energy E and the density [Eq. (6)] of the system become exactly equal to the energy and the density of non-interacting spinless fermions, respectively. For our showcase of few-bosons systems (N = 2 to N = 5) in a harmonic trap with frequency one, We start our investigation with the one-body density as a function of the interaction strength λ [ Fig. 1 For comparatively weak repulsion, the density is clustered at the center of the trap, but becomes flatter and broader when λ increases. For stronger repulsion, the density gradually acquires modulations and the number of humps finally saturates to the number of bosons in the system; four humps for N = 4 bosons are clearly visible when the interaction strength goes above λ ~ 10. The emergence of N maxima in the density indicates that the TG regime is approached. The density modulations/humps are more pronounced in the center of the trap, where the potential is close to zero. For larger distances from the origin, the humps in the density are less pronounced due to the non-zero value of the confinement potential. Importantly, the outermost density modulation also becomes less pronounced if the number of particles N is increased. See also Appendix B for a direct comparison of the relative height of innermost and outermost peaks for different particle numbers N.
We note that the density's maxima in the Tonks-Girardeau regime are distinct but not isolated. We also observe that, once the TG regime is reached, the density does not broaden further with increasing values of λ for all particle numbers. We also provide a direct comparison with the ground state properties of non-interacting fermions computed with the multiconfigurational time-dependent Hartree method for fermions (MCTDHF), see Fig. 2. We note that the results for non-interacting fermions can be obtained analytically, i.e., here, we use the heavy MCTDHF method only for the sake of computational convenience.
We now move to discuss the two-body densities ρ (2) For stronger repulsion a so-called "correlation hole" in the two-body density forms on the diagonal, ρ (2) Fig. 3(a) for λ = 10 and λ = 30]. The probability of finding two bosons at the same position tends towards zero. In the limit of infinite repulsion the correlation hole persists in ρ (2) . In analogy, however, to the boundedness of the energy as a function of the interaction strength, the width of two-body density on its anti-diagonal [ρ (2) (x, −x)] is also bounded, i.e., the spread of ρ (2) converges in the fermionization limit when λ → ∞.
Similar to the one-body density, the maxima which are formed in the off-diagonal of the two-body density are distinct but not isolated [see Fig. 3(a) for λ = 10 and λ = 30].
We infer that the correlation hole along the diagonal and the confined spread are the unique signatures of the two-body density of a fermionized state.
crystallization. For bosons with dipole-dipole interactions, crystallization occurs when the long-range tail of the interaction [see Eq. (4)] becomes dominant 34 : the bosons form a lattice structure which allows them to minimize their mutual overlap. To characterize crystallization we analyze the one-body and two-body density for bosons with dipolar interaction of strength g d . We choose the cut-off parameter α = 0.05 in Eq. (4) such that the effective interaction features the same physical beahavior as the "real" dipolar interaction that additionally contains a contact-interaction term (see Ref. 34 for a direct comparison).
We plot the one-body density of N = 4 bosons as a function of g d in Fig. 1(c,d). The system is condensed at the center of the trap for small g d . As g d increases, the density starts to exhibit a four-hump structure (see Fig. 1(c,d) for ∈ ∼ ∼ g [ 1,5] d ) similar to the density observed for the fermionization of bosons with contact interactions [ Fig. 1(a)]. www.nature.com/scientificreports www.nature.com/scientificreports/ This attempted fermionization results from a dominant contribution of the short-range part of the dipolar interaction potential, see also Ref. 18 . However, this fermionization-like behavior is only a precursor to the crystal transition that takes place when the long-range nature of the interaction starts to dominate the physics of the system for larger interaction strengths [ Fig. 1(c,d) for g 10 d  ]. For crystallized dipolar bosons at sufficiently large g d , the value of the density at its minima between the humps tends to zero while the spreading of the density profile diverges as g d increases, see Fig. 4. At g d = 30.0, we observe four well-isolated peaks heralding the crystallization of the N = 4 bosons. We collect results for other numbers of bosons (N = 2, 3, 5, 6) with dipole-dipole interactions -including the relative height of the peaks in the density that shows that the peaks are well-isolated in comparison to particles with contact interactions -in Appendix B. A comparison of MCTDHB results with exact diagonalization is shown in Appendix A.
We now analyze the two-body density for dipole-dipole interactions [ Fig. 3(b)]. For small interaction strength, g d = 0.1, the atoms are clustered together at the center of the trap. As g d increases, a correlation hole develops: ρ (2) (x, x) tends to zero [ Fig. 3(b) for g d ≥ 1]. Thus, due to the long-range interaction, the probability of finding two bosons in the same place is strongly reduced. In the crystalline phase [ Fig. 3(b) for g d ≥ 10]: the bosons escape their spatial overlap entirely and even the off-diagonal peaks of ρ (2) become isolated. We term this behavior the formation of an off-diagonal correlation hole. For crystallized bosons, the spread of the anti-diagonal of the two-body density, ρ (2) (x, −x), is diverging as g d is increasing [compare Fig. 3(b) for g d = 10 to Fig. 3(b) for g d = 30].
We assert that the correlation hole along the diagonal and the off-diagonal and the unbounded spreading are the unique signatures of the two-body density of a crystalline state of dipolar bosons.
Sorting crystallization from fermionization. We now discuss how to distinguish fermionized from crystallized many-body states. One clear distinction is given by the spread of the one-and two-body densities: for bosons with contact interactions it is bounded, while for bosons with dipole-dipole interactions it diverges as a function of the interaction strength. We assert that, (1) the bounded spreading of the density for contact www.nature.com/scientificreports www.nature.com/scientificreports/ interactions is a consequence of the bounded energy as the interaction strength tends to infinity. Similarly, we assert, (2) that the unbounded spreading of the density for dipole-dipole interactions is a consequence of the unbounded energy as the interaction strength g d tends to infinity. To validate the assertions (1) & (2), we quantify the spreading of the density as a function of interactions and plot the position of its outermost peak as a function of the interaction strength in Fig. 5(a) and compare it to the energy in Fig. 5 From fitting the energy in Fig. 5(b) we can infer that the energy as a function of contact interaction strength approaches the fermionization limit exponentially, a power law does not fit as accurately the data. For very large interactions, in the limit of λ −1 → 0, our results are in agreement with the analysis in Ref. 67 , see Appendix B. For dipolar interactions, the growth of the energy as a function of the interaction strength is fitting well to a power law.
We thus conclude that the crystalline phase can be distinguished from the TG regime gas by virtue of the behavior of its density profile as a function of the strength of the interparticle interactions: The width of the density distribution converges for an increasing strength of contact interactions, but it continuously spreads for an increasing strength of long-range interactions [compare Fig. 1(b) with Fig. 1(d) as well as Fig. 5(a) with Fig. 5(b)]. In Appendix B, we demonstrate that the exponent of the power law of the spreading of the density as a function of the strength of the interaction is independent of the particle number N.
In the case of long-ranged interactions, the unbounded spreading of densities as a function of increasing interaction strength and the formation of well-isolated peaks are in sharp contrast to the bounded spreading of densities and the non-isolated peaks in the case of contact interactions in the TG regime [cf. Figs 1 and 5(a-b)].
We now turn to analyze the eigenvalues of the reduced one-body density matrix, the so-called natural occupations 60 , as a function of the interaction strength between the particles [ Fig. 5(c,d)]. As expected 8,33,34 , when the value of the interaction strength increases, the occupation of the first natural orbital decreases while the other orbitals start to be occupied. For contact interactions, mostly one natural occupation, n 1 , dominates, while the other occupations n k , k > 1 remain comparatively small even for large values of λ: depletion emerges as the fermionized state is reached [ Fig. 5(c)], see also Ref. 8 . For long-range interactions, however, all occupations ρ k (NO) for k ≤ N contribute on an equal footing for large values of g d . This full-blown N-fold fragmentation emerges as the crystal state is reached [ Fig. 5(d)], see also Ref. 34 . In the crystal state the bosons behave similar to distinguishable particles 68 , and the particle statistics does not influence the physical observables considered in Fig. 5(b,d): the www.nature.com/scientificreports www.nature.com/scientificreports/ energy and the natural occupations for bosons and fermions converge to the same values. Thus, the finding of Ref. 68 for two particles may be extended to larger number of particles.
The emergence of complete fragmentation is a consequence of long-ranged interactions and in sharp contrast to the emergent depletion in the case of contact interactions.

conclusions
In this paper we highlight the key characteristics of the many-body wavefunction that reveal the difference between the fermionized bosons with contact interactions and crystallized bosons with dipolar interactions.
In the case of fermionization, the one-(two-)body density shows a modulation with a number of maxima corresponding to the number of particles. The maxima are confined but not completely separated. The incomplete separation is a consequence of the representability of momentum distribution of fermionized bosons using a basis set: infinitely many basis states are necessary to accurately resolve the cusp -a fact that is reflected by the depletion of the state which we quantified by the eigenvalues of the reduced one-body density matrix. We found that the peaks in the density as well as the energy as a function of the interaction strength approach the fermionization limit exponentially.
In the case of crystallization, the one-(two-)body density shows well-separated peaks whose distances diverge as a function of the interaction strength as a power law. This completed separation is the consequence of the formation of a Mott-insulator-alike many-body state where the "lattice potential" is replaced by the long-ranged interparticle interactions and the "lattice constant" is dictated by the strength of the interparticle interactions.
We close by stating that all the signatures that distinguish crystalline bosons from fermionized bosons can be measured experimentally using single-shot absorption imaging [69][70][71][72][73] . From experimental absorption images, the one-body and two-body density are available as averages of many single-shot images. Thus, a direct verification of our results for the spread of the one-body and two-body density can be performed. Furthermore, Refs. 34,63 suggest that the natural occupations can be inferred from the integrated variance of single-shot images, at least at zero  Fig. 2 for the fitting parameters A N and B N . Furthermore, we assess the convergence of the spread of the density as a function of the interaction strength to the spread of the density of the non-interacting fermionic system, see arrows labeled "2F","3F","4F", and "5F" in Fig. 2.
To compare our results for the energy in the fermionization limit to analytical predictions for very large contact interaction strengths in Ref. 67 , we plot the energies as a function of −λ −1 in Fig. 7. We find that our results are consistent with the linear limit for the energy as a function of −λ −1 of Ref. 67 .
We now turn to the relative height of the innermost and outermost peak(s),    bosons with contact interactions (lines, top, to bottom, respectively). The relative peak height is consistently smaller for the outermost peak as compared to the innermost peak in the density for all interaction strengths depicted. All quantities shown are dimensionless. www.nature.com/scientificreports www.nature.com/scientificreports/ in the density. Here, ρ max refers to the value of the density ρ(x) at the peak position and ρ min refers to the value of the density ρ(x) at the position of the minimum to the left to the considered peak. See Fig. 8 for a plot of Δρ(x) for N = 2, 3, 4, 5 bosons. It is clearly seen that, for fixed N, the outermost peaks' relative height is much smaller than the relative height of the innermost peaks.

Dipolar interactions.
Here, we assess the validity of the power-law-like unbounded spreading of the density as a function of the strength of dipole-dipole interactions, that we have shown in Fig. 4 of the main text for N = 4 particles. In Fig. 4 We plot the spread of the density for N = 2, 3, 4, 5, 6 dipolar bosons obtained with MCTDHB with M = 16, 16, 22, 28, 26 orbitals, respectively, and fit it with a power law C g N x D N . We find that the exponent in the power law is almost identical for all particle numbers studied here, i.e., D N ≈ 0.16 for N = 2, 3, 4, 5, 6.
We now discuss the relative peak height Δρ(x), see Eq. (B1), of the outermost peak as a function of the dipolar interaction strength, see Fig. 9 for a plot for N = 2, 3, 4, 5, 6. As hinted by Fig. 1(c) in the main text, the relative peak height for the case of the dipole-dipole interactions converges towards unity as the strength of interactions g d increases, because the values of the minimum, ρ min in Eq. (B1), tends to zero: the peaks in the crystal state are well-isolated in comparison to the peaks in the fermionization limit for bosons with contact interactions (compare magnitude of relative peak heights in Figs. 8 and 9).