S Bera

B Chakrabarti

A Gammal

M C Tsatsos

M L Lekala

B Chatterjee

C Lévêque

Axel U J U J Lode

Sorting Fermionization from Crystallization in Many-Boson Wavefunctions

Keywords: Schrödinger equation, fermionization, crystallization, dipolar interactions, many-body physics

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The physics of the ultracold Bose gas in one spatial dimension is strongly different from that of its threedimensional counterpart [1,2]. In one spatial dimension, quantum effects are much more prominent. When the interactions are strong, quantum fluctuation are enhanced. Experimentally, in cold atom systems, the dimensionality can be manipulated using a tight transversal confinement that essentially freezes the radial degrees of freedom [3,4]. Such quasi-one dimensional systems display intriguing physics that cannot be realized for threedimensional systems: Fermionization, occurs for strongly interacting bosons with contact interactions [5][6][7][8][9][10][11] and crystallization emerges for sufficiently strongly interacting bosons with dipole-dipole interactions [12][13][14][15][16][17]. For bosons with contact interactions, fermionization leads to the formation of the Tonks-Girardeau (TG) gas when the interaction strength tends to infinity. This is a consequence of the Bose-Fermi mapping [7,[18][19][20] which implies that strongly interacting bosons and non-interacting spinless fermions have the same one-body density in position space. With increasing interaction strength, not only the density, but also the energy of fermionized bosons saturates to the energy of non-interacting fermions.

In the case of dipolar interactions, the remarkable phenomenon of crystallization occurs when the interaction strength is moderately large. Bosons interacting via a dipole-dipole interaction potential have become the primary cold atom system to investigate the many-body physics triggered as consequence of long-range interactions [21][22][23][24][25][26]. The long-ranged and anisotropic nature of the dipolar interaction potential results in a variety of interesting effects and phenomena [27], like crystallization in one-dimensional systems [12][13][14][15][16][17], that are completely different from the emergent phenomena in the case of strong contact interactions. Crystallization is a consequence of the repulsive and long-ranged tail of the dipolar interactions dominating the physics [28,29]: the bosons maximally separate and minimize their spatial overlap. Unlike in the fermionization of bosons with strong contact interactions, the energy of crystallized bosons does not saturate. We note here that it is formally possible to define and measure an order parameter which is a function of the eigenvalues of the reduced one-body density matrix that allows to unequivocally identify the crystal phase of one-dimensional dipolar bosons [29]. Further-more, we note that the formation of a crystal state is a generic feature of many-body systems of particles with long-ranged interactions. Fermions with long-ranged interactions, for instance, form a so-called Wigner crystal [30].

In this work, we analyze the differences between fermionized and crystallized bosons' wavefunctions using the energy as well as the one-body and two-body reduced density matrix. We demonstrate how fermionization can be distinguished from crystallization by quantifying the (experimentally accessible) spread of the onebody and two-body densities. The different spreading characteristics of the one-body and two-body densities for fermionized as compared to crystallized bosons are a direct consequence of the different behavior of the energy as a function of interaction strengths. For dipolar interactions the energy as a function of interaction strength is unbounded; this is in stark contrast to the bounded energy as a function of interactions for contact interactions. Current experimental setups, for instance, for Erbium [31], Erbium-Erbium molecules [32], or Sodium-Potassium molecules [33], enable the experimental exploration of systems with dominant dipole-dipole interactions needed for probing the physics of crystallization.

The (momentum) densities of bosons with dipolar interactions have been compared to those of fermions with dipolar interactions in Ref. [17]. Going beyond Ref. [17], we compare and quantify the spreading of the full density matrices of bosons with dipolar interactions to the spreading of the density matrices of bosons with contact interactions. We note that Ref. [34] discusses and compares the physics of spin-1/2 fermions with contact and with long-ranged interactions in lattices using a Hubbard description. Our work complements the findings in Ref. [34] by providing a comparison of single-component ("spin-0") bosons with contact and long-ranged dipolar interactions.

Fermionization and crystallization entail the breakdown of mean-field approaches like the time-dependent Gross-Pitaevskii equation [START_REF] Pethick | Bose-Einstein Condensation in Dilute Gases[END_REF][START_REF] Bogoliubov | Selected Works II: Quantum and Statistical Mechanics[END_REF][START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF]. We thus use the multiconfigurational time-dependent Hartree for bosons (MCTDHB) [START_REF] Ofir | [END_REF]39] and fermions (MCTDHF) [40] methods implemented in the MCTDH-X software package [40][START_REF] Lode | MCTDH-X:The time-dependent multiconfigurational Hartree for indistinguishable particles software[END_REF][START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF][43] to compute the ground state of the few-particle Schrödinger equation.

We illustrate our findings with computations for N = 4 bosons in a parabolic trapping potential and trace the complete range of dipolar and contact interaction strengths by obtaining highly accurate results with MCT-DHB.

This paper is structured as follows: in Sec. II, we discuss the Hamiltonian and quantities of interest, in Sec. III, we introduce the numerical method, MCT-DHB, that we use for obtaining solutions of the few-body Schrödinger equation, in Sec. IV we analyze fermionized and crystallized bosonic few-body states and discuss how they can be sorted from each other and we conclude our paper in Sec. V. Results for other observables and an as-sessment of the accuracy of our computations with the exact diagonalization and MCTDHF approaches are collected in the Appendices A and B.

II. HAMILTONIAN, ONE-AND TWO-BODY DENSITY

In order to discuss crystalline and fermionized bosons, we consider the time-dependent many-body Schrödinger equation,

i∂ t |Ψ = Ĥ|Ψ . (1)
Here, |Ψ is the many-body state and the N -particle Hamiltonian Ĥ in dimensionless units [START_REF]considered bosons and L a conveniently chosen length scale[END_REF] reads,

H = N i=1 - 1 2
∂ 2 ∂x 2 i + V (x i ) + i<j W (x i -x j ), (2)
where we set V (x i ) = 1 2 x 2 i to be the external harmonic trap and W (x i -x j ) the interaction potential. All quantities are dimensionless and expressed in harmonic oscillator units. For contact interactions

W (x i -x j ) = λδ(x i -x j ), (3)
where λ is the interaction strength determined by the scattering length a s and the transverse confinement frequencies [START_REF] Olshanii | [END_REF]. For long-ranged dipolar interactions

W (x i -x j ) = g d |x i -x j | 3 + α , (4)
where g d is the dipolar interaction strength and α is a short-range cut-off to avoid the divergence at x i = x j . We restrict ourselves to N = 4 interacting bosons for all our calculations and consider repulsive interactions, λ > 0 and g d > 0, exclusively. Our simple approximation to the one-dimensional dipole-dipole interaction potential in Eq. ( 4) is justified for the moderate to large interaction strengths 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 [46]. 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 dipoledipole interaction augmented with an additional contact interaction potential, see Ref. [29]. A rigorous discussion of the dipole-dipole interaction potential in one and two spatial dimensions can be found in Ref. [46]; here, for the sake of computational complexity, we will focus on quasi-one-dimensional systems.

To investigate the stationary properties of the system, we propagate the wavefunction in imaginary time using the MCTDH-X software [40][START_REF] Lode | MCTDH-X:The time-dependent multiconfigurational Hartree for indistinguishable particles software[END_REF][START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF][43] to solve the MCTDHB equations of motion [START_REF] Ofir | [END_REF]39]; thereby the system relaxes to the ground state.

In the following we discuss the reduced one-body density matrix, defined as

ρ (1) (x, x ) = Ψ| ψ † (x ) ψ(x)|Ψ . (5)
Its diagonal,

ρ(x) ≡ ρ (1) (x, x = x) (6)
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 ρ

(N O) i
of the reduced onebody density matrix ρ (1) in Eq. (5). For this purpose, we write ρ (1) in its eigenbasis:

ρ (1) (x, x ) = i ρ (N O) i Φ * i (x )Φ i (x). ( 7
)
The eigenvalues ρ

(N O) i
and eigenfunctions Φ i (x) are referred to as natural occupations and natural orbitals, respectively. If only a single eigenvalue ρ

(N O) i
is macroscopic, then the state |Ψ describes a Bose-Einstein condensate [47]. The case when multiple eigenvalues ρ

(N O) i
are comparable to the number of particles N is referred to as fragmentation [START_REF] Ofir | [END_REF]39,[START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF]43,[48][49][50].

In the following, we will also use the two-body density ρ (2) to characterize crystallization and fermionization. It is defined as

ρ (2) (x 1 , x 2 ) = Ψ| ψ † (x 1 ) ψ † (x 2 ) ψ(x 1 ) ψ(x 2 )|Ψ . (8)
The two-body density quantifies the probability to detect two particles at positions x 1 and x 2 .

III. NUMERICAL METHOD

The computation of the exact many-body wave function is a difficult problem. To attack the time-dependent many-body Schrödinger equation, Eq. ( 1), we expand the many-body wavefunction |Ψ of N interacting bosons in a complete set of time-dependent permanents | n; t = |n 1 , ..., n M ; t with at most M single-particle states or orbitals. The MCTDHB ansatz for the many-body wave function is thus

|Ψ(t) = n C n (t)| n; t . (9)
Here, the permanents | n; t are symmetrized bosonic many-body states that are also referred to as "configurations". The sum in Eq. ( 9) runs on all configurations n of N particles in M orbitals. The number of perma-

nents and coefficients C n (t) is N conf = N + M -1 N .
In the second quantized representation the permanents are given as

|n; t = |n 1 , ...n M ; t = M i=1   b † i (t) ni √ n i !   |vac . ( 10
)
Here b † k (t) is the bosonic creation operator which creates a boson in the time-dependent single particle state φ k ( r, t). Eq. ( 9) 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 singleparticle states M . A set of coupled equations of motion for, both, the time-dependent expansion coefficients C n (t) and the time-dependent orbitals φ k ( r, t) are obtained by requiring the stationarity of the action of the time-dependent Schrödinger equation [START_REF] Ofir | [END_REF]39] under variations of C n (t) and φ k ( r, t). Using MCTDHB, both, the coefficients and orbitals are variationally optimized [START_REF]Geometry of the time-dependent variational principle[END_REF]. 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 Ref. [START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF][START_REF] Lode | [END_REF] 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). We solve the set of coupled MCTDHB equations simultaneously for the ground state using imaginary time propagation [40][START_REF] Lode | MCTDH-X:The time-dependent multiconfigurational Hartree for indistinguishable particles software[END_REF][START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF][43].

IV. 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. IV A and Sec. IV B, respectively. Thereafter, we investigate how to sort the one, fermionization, from the other, crystallization, in Sec. IV C. 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.

A. 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 system N = 2 to N = 5 bosons in a harmonic trap with frequency one, V (x) = 1 2 x 2 , the limiting value is thus E N λ→∞ = N 2 2 . We start our investigation with the one-body density as a function of interaction strength λ [Fig. 1(a-b)]. 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 provides a direct comparison with the ground state properties of non-interacting fermions computed with the multiconfigurational timedependent Hartree method for fermions (MCTDHF), see Fig. 2.

We now move to discuss the two-body densities ρ (2) of bosons with contact interactions [see Fig. 3(a)]. For weak interaction strength, λ = 0.1, the bosons are clustered near the center, i.e, x 1 = x 2 = 0 [Fig. 3(a)]. As the interaction strength increases, ρ (2) spreads out to the off-diagonal (x 1 = x 2 ) while the diagonal (x 1 ∼ x 2 ) is depleted [see Fig. 3(a) for λ = 1].

For stronger repulsion a so-called "correlation hole" in the two-body density forms on the diagonal, ρ (2) (x, x) → 0 [see 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 antidiagonal [ρ (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.

B. Crystallization

For bosons with dipole-dipole interactions, crystallization occurs when the long-range tail of the interaction [see Eq. ( 4)] becomes dominant [29]: 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 V ef f = 10.0 -10.0 g d x 3 +α dx = 10.0 -10.0 δ(x)dx = 1.0. We plot the one-body density of N = 4 bosons as a function of g d in Fig. 1(c) and (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) and (d) for g d ∈ [∼ 1, ∼ 5]) similar to the density observed for the fermionization of bosons with contact interactions [Fig. 1(a)].

This attempted fermionization results from a dominant contribution of the short-range part of the dipolar interaction potential, see also Ref. [17]. 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) and (d) for g d 10]. For crystallized dipolar bosons at 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 offdiagonal 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 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.

C. 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 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(b), for N = 4.

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. For very large interactions, in the limit of λ -1 → 0, our results are in agreement with the analysis in Ref. [53], 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.

Indeed, the comparison of Fig. 5(a) and Fig. 5(b) corroborates our assertions 1) & 2), and holds for different number of particles.

We thus conclude that the crystalline phase can be distinguished from the Tonks 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. Fig. 1 and Fig. 5(a)-(b)].

We now turn to analyze the eigenvalues of the reduced one-body density matrix, the so-called natural occupations [47], as a function of the interaction strength between the particles [Fig. 5(c)-(d)]. As expected [8,28,29], 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)], cf. also Ref. [8]. For long-range interactions, however, all occupations ρ

(N O) k
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. [29].

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.

V. 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 [54][55][56][57][58]. 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. [29,50] suggest that the natural occupations can be inferred from the integrated variance of single-shot images, at least at zero temperature. It is, of course, an open question how thermal fluctuations affect the variance in absorption images and up to which temperature it is still possible to determine the fragmentation of the system.

under the DST Inspire Faculty fellowship. AUJL and CL acknowledge financial support by the Austrian Science Foundation (FWF) under grant Nos. P 32033 and M 2653, respectively, and the Wiener Wissenschafts-und TechnologieFonds (WWTF) project No MA16-066 . The results of the manuscript have been obtained with MCTDHB with M = 12, 14, 20, 22 orbitals for N = 2, 3, 4, 5 bosons, respectively, with a contact interaction strength up to λ = 1000. In Fig. 2 of the main text, our results are consistent with fits of an exponential function A N (exp(-λ/B N ) -1), see caption of 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. [53], 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. [53].

We now turn to the relative height of the innermost and outermost peak(s), ∆ρ(x) = ρ max -ρ min ρ max + ρ min (B1) 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 N g D N d . 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. 1c) 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 and9).

  (b) for g d = 10 to Fig. 3(b) for g d = 30].

FIG. 1 . 5 FIG. 2 . 6 FIG. 4 .FIG. 5 .

 152645 FIG. 1. One-body density of four bosons as a function of contact [(a),(b)] and dipolar [(c),(d)] interparticle interaction strength. For contact interactions, the density becomes flatter and broader as the repulsion increases [panel (a) and (b) for λ ≤ 1]. For even larger interaction strengths [panel (a) and (b) for λ 10], four distinct but not isolated peaks appear and the density gradually converges to the density of four non-interacting fermions as λ → ∞. Due to this convergence, the spread of the density seizes to increase [panel (d)]. For dipolar interactions, the one-body density is clustered at the center of the trap for small interactions [panels (c),(d) for g d 1]. As g d increases, the density develops a fourfold splitting [panel (c) and (d) for g d 1]. As a function of increasing interaction strength, the spread of the density continues to increase [panel (d)] and the fourfold spatial splitting intensifies to form four completely isolated peaks in the density for sufficiently strong dipolar interactions: crystallization emerges [panels (c),(d) for g d 10]. All quantities shown are dimensionless.

FIG. 6 .

 6 FIG. 6. Comparison of MCTDHB and ED for N = 4 bosons with dipole-dipole interaction strength g d = 30. The plot shows the relative error in energy with respect to an MCT-DHB computation with M = 32 orbitals as a function of the number of orbitals for the ED and MCTDHB approaches. Due to the variationally optimized basis in MCTDHB computations it features a much smaller error for any number of orbitals. All quantities shown are dimensionless.

5 FIG. 7 . 5 FIG. 8 . 6 FIG. 9 .

 575869 FIG.7. Energy as a function of inverse interaction strength, -1/λ, for N = 2, 3, 4, 5 (bottom to top curve, respectively) bosons. Our results are consistent with the analysis in Ref.[53]: the energy linearly converges to the fermionization limit, i.e., when -λ -1 → 0. All quantities shown are dimensionless.

ACKNOWLEDGMENTS S. Bera acknowledges DST (Govt. of India) for the financial support through INSPIRE fellowship [2015/IF150245] and Rhombik Roy for some useful discussions. B. Chakrabarti acknowledges FAPESP (grant No. 2016/19622-0). AG and MCT acknowledge FAPESP and AG thanks CNPq for financial support. B. Chatterjee acknowledges financial support from the Department of Science and Technology, Government of India

Appendix A: Comparison of MCTDHB and exact diagonalization

In this Appendix, we demonstrate that MCTDHB yields solutions to the Schrödinger equation at a larger accuracy as compared to the exact diagonalization approach (ED). As is conventional, we use the eigenfunctions of the non-interacting system as the single-particle basis states for the ED. We solve the same system as shown in Fig. 1c) and d) for an interaction strength of g d = 30 and compare the energies obtained with MCT-DHB and ED, see Fig. 6. Due to the variationally optimized single-particle basis in MCTDHB computations it features a much smaller error than the ED computations with an unoptimized single-particle basis for the same number of orbitals. This observation is in agreement with other works that benchmark the MCTDHB and the MCTDHF approaches against ED, see Ref. [START_REF] Lode | [END_REF] and Ref. [40], respectively.

Appendix B: Different particle numbers

In this Appendix, we corroborate our results in the main text by studying different particle numbers.


Budhaditya Chatterjee

Jörg Schmiedmayer
email: schmiedmayer@atomchip.org

Camille Lévêque
email: camille.leveque@tuwien.ac.at§axel.lode@univie.ac.at

Axel U J Lode

Unveiling Emergent Crystal Orders of Incommensurate Dipolar Bosons in One-Dimensional Lattices using Full Distribution Functions

We explore the ground-states of a few dipolar bosons in optical lattices with incommensurate filling. The competition of kinetic, potential, and interaction energies leads to the emergence of a variety of crystal state orders with characteristic one-and two-body densities. We probe their transitions and construct the emergent state diagram as a function of the dipolar interaction strength and the lattice depth. We demonstrate that the crystal state orders can be observed using the full distribution functions of the particle number which we extract from simulated single-shot images.

The experimental realization of stable Bose-Einstein condensates (BECs) of dipolar atoms [1][2][3][4][5][6] and molecules [5,6] provides new perspectives to study the phase transitions in correlated quantum systems [7,8]. The dipole-dipole interactions are of anisotropic and long-range nature causing a plethora of new phenomena absent in conventional BEC, e.g., directional elongation [9][10][11] and geometric stabilization [9,10,[12][13][14]. The lower dimensionality results in additional physical features: p-wave superfluidity [15,16], Luttinger-liquidlike behavior [17][18][19][20] and anisotropy in curved geometries [21][22][23].

Trapped atoms in optical lattices serve as quantum simulators for condensed matter systems [24][25][26][27][28][29][30][31][32][33][34], exploiting the precise experimental tunability of these quantum systems. Even few-particle systems have been experimentally studied providing a bottom-up access to manybody physics [35][36][37]. The interplay between anisotropic long-range interaction and contact interaction in dipolar bosons results in the emergence of new phases: Besides the usual superfluid and Mott-Insulator phases, a density wave phase (DW ) [38,39], characterized by an alternate filling of lattice sites and a supersolid phase [40][41][42][43][44], with coexistent DW structure and superfluidity, and more exotic phases such as a Haldane insulating phases [39,45], checkerboard phases [38,46] and Mott solids [47] were predicted. Dipolar atoms with strong interactions feature a crystal state (CS) [17, 21, 22, 48-50, 52, 53].

In this letter, we investigate the plethora of crystal orders that emerge in one-dimensional lattices incommensurately filled with dipolar bosons and put forward protocols to detect them experimentally. Crystal states with different orders emerge for different strengths of the dipolar interaction, due toa competition between kinetic, potential, and dipolar interaction energies. We demonstrate that all the emergent crystal orderings can be identified using the full distribution functions of the positiondependent particle number operator extracted from absorption or single-shot images. We study a few dipolar atoms in a one-dimensional optical lattice and determine the ground-state properties numerically with MCTDHB [59,60] implemented in the MCTDH-X software [61][62][START_REF] Lode | MCTDH-X:The time-dependent multiconfigurational Hartree for indistinguishable particles software[END_REF]. Arbitrarily strong interacting systems can be investigated with MCTDHB [START_REF] Fischer | [END_REF][65][66][67][68][69][70][71] allowing us to explore the full range of interaction strengths and observe two distinct crystal orders. We specify two observables to characterize these distinct orders of the crystal state and compute the state diagram, i.e., the finite-size analog of the phase diagram in the thermodynamic limit. We establish experimental protocols for the detection of the state diagram based on the full distribution function, the variance, and a binning analysis of single-shot images [52,[72][73][74].

We consider polarized, dipolar bosons in a quasi-onedimensional lattice potential, V ol = V sin 2 (κx), with a depth V and a wave-vector κ. A tight transversal confinement with a characteristic length a ⊥ prevents an excitation into the transverse direction and ensures the quasi-one-dimensionality of the system. The interaction is purely dipolar,

V int (x i -x j ) = g d |xi-xj | 3 +α ; the dipo- lar interaction strength, g d , is g d = d 2
m /4π 0 for electric dipoles and g d = d 2 m µ 0 /4π for magnetic dipoles, d m is the dipole moment, 0 the vacuum permittivity, and µ 0 the vacuum permeability [75]. For large separations,

|x i -x j |
a ⊥ , we recover the far-field dipole-dipole interaction ∼ 1/r 3 and for small separations, |x i -x j | a ⊥ , the transversal confinement introduces a short-scale cutoff α ≈ a ⊥ 3 [49,76,77]. We make all quantities dimensionless by expressing them in terms of the lattice recoil energy for particles of mass M , E R = 2 κ 2 /2M . In this work, we set the cutoff α = 0.05 which corresponds to an aspect ratio γ = 25.6 and consider N = 8 bosons. We impose hard-wall boundary conditions to restrict our lattice to S = 5 sites. Our setup can be generalized for larger ensembles.

We analyze possible orderings in the crystal state with the density ρ(x) = Ψ| Ψ † (x) Ψ(x)|Ψ as a function of the interaction strength g d [Fig. 1(a)]. For g d = 0, we obtain a pure superfluid state (SF ). As g d increases, the incommensurate setup implies the absence of a pure Mott-Insulator state (M I) since; a SF coexists with the M I [78]. This coexistence of M I and SF (henceforth SM I) is seen in the structure of ρ(x): the two outer wells have a smaller population than the central wells; N = 5 atoms in a M I state coexist with N = 3 bosons in the SF fraction. The SF localizes in the central wells because this minimizes the kinetic energy in our setup with hard wall boundaries.

The crystal transition occurs at g d ≈ 1: due to repulsive dipolar interactions, the bosons avoid each other and minimize their overlap to save interaction energy. The one-body density [Fig. 1(a)] exhibits a splitting of the density in the doubly occupied central wells -a signature of onset of the crystal transition. In terms of bosons occupation per site,the crystal state is ordered as 1, [11], [11], [11], 1. Here, a two-hump density in doublyoccupied lattice sites in the crystal state is represented as [11]. Since the double occupation of the three central wells is a result of the kinetic energy term in the Hamiltonian, we term this state "kinetic crystal state" (KCS).

Upon a further increase of the interaction strength g d , the repulsive tail of the long-range dipolar interactions overcomes the kinetic energy and makes the double occupation of adjacent sites energetically unfavorable [Fig. 1(a)]: instead of nearest neighbors, next-nearest neighbors are doubly occupied in a density-wave-ordered structure at g d 3. For our setup, a pure density-waveorder would have the occupations 2, 1, 2, 1, 2. Due to the strong dipolar interaction, the bosons density in doubly occupied sites is spatially split and a density-wave crystal state (DW CS) with occupations [11], 1, [11], 1, [11] is obtained. The KCS and the DW CS are among possible crystal orders and hence subsets of a general crystal state CS. Here, and henceforth, we use the label CS for a crystal state in either the KCS or the DW CS arrangement.

For very large values of the interaction g d > 15, the long-range interactions overwhelm the lattice potential [SI [88], Sec. S2]. Consequently, the DW CS gradually transmutes to an equispaced crystal state where the bosons are not localized around the potential minima anymore (not shown). See SI [88], Sec. S3 for other possible crystal orders at different particle numbers and lattice sizes.

The various crystal arrangements of the atoms with respect to the lattice can be obtained from a purely classical model (SI [88], Sec. S3). However, the many-body calculations we present are necessary to capture the quantum properties of the crystal state: many modes do contribute to the quantum field of the crystal state (cf. Fig. 1) and thus its properties can be assessed only through the analysis of a realistic model of the many-body wavefunction. Moreover, the classical description does not contain a SF or M I state so the transition from a non-crystal to a crystal order cannot be obtained from a classical model.

To unveil the mechanisms underlying the localization and coherence in the emergent states, we now discuss the two-body density ρ2 (x 1 , x 2 ) = Ψ| Ψ † (x 1 ) Ψ † (x 2 ) Ψ(x 1 ) Ψ(x 2 )|Ψ for characteristic values of g d (Fig. 2). Unlike the particle arrangement corresponding to the one-body density [Fig. 1(a)], our analysis of the two-body density exhibits true quantum manybody properties of the crystal state: the two-body densities in Fig. 2(b)-(d) -unlike for classical, semi-classical, and mean-field models -cannot be represented as a product of one-body densities.

For small interaction strength in the SF (g d = 0.0005) the maxima of ρ (2) are nearly uniformly distributed where both arguments x 1 and x 2 are in the vicinity of a minimum of the lattice potential [Fig. 2(a)]. For a larger interaction, g d ≈ 0.1, a partial depletion along the diagonal of ρ (2) (x 1 ≈ x 2 ) occurs. This depletion results from the formation of an M I that coexists with a SF fraction [Fig. 2(b)]. At even stronger interactions, g d = 1.5, the bosons in the doubly-occupied central wells crystallize forming the KCS. Here the diagonal of ρ (2) is completely depleted, ρ (2) ≈ 0 for x 1 ≈ x 2 : a correlation hole is formed, the probability of detecting two bosons at same position (x 1 = x 2 ) becomes negligible. The split maxima in the three central wells result from the on-site interaction-driven splitting 2 → [11]. For g d ≈ 6, the system is in the DW CS [Fig. 2(d)]. A split inter-site structure of the two-body density ρ (2) is now present for every odd site of the lattice potential. The diagonal depletion is larger compared to the KCS [compare Fig. 2(c) and (d) for x 1 ≈ x 2 ]. The non-uniform distribution of the maxima of ρ (2) clearly signals the DW CS.

When the interactions increase even further, the DW CS -to minimize the interaction energy -transforms into a crystal state (not shown) where the distribution of the maxima in both ρ(x) and ρ (2) (x 1 , x 2 ) becomes equispaced and independent of the underlying lattice potential. This crystal state for very strong interactions (g d 15) is different from the SM I, the KCS, and the DW CS, since the distribution of maxima of the densities is not dictated by the minima of the lattice potential.

The transition from the SM I to the KCS and the DW CS depends also on the lattice depth V . To address the role of the lattice depth V and the interaction strength g d , we construct the state diagram of the system as a function of g d and V . The SM I → CS transition can be determined from a state order-parameter constructed from eigenvalues of the one-body reduced density matrix (RDM) as defined in Ref. [52]. The first order RDM is defined as

ρ(1) (x, x) = Ψ| Ψ † (x) Ψ(x)|Ψ = i λ i ϕ * i (x)ϕ i (x ). ( 1
)
By diagonalizing ρ(1) , the eigenvalues λ i (natural occupations) and eigenfunctions ϕ i (x) (natural orbitals) are obtained. The values λ i determine the degree of condensation or fragmentation of a state. If a single natural occupation is macroscopic (λ 1 ≈ N ), the system is condensed [79] and if multiple λ i s are macroscopic, the system is fragmented [80,81].

The natural occupations are shown as a function of the interaction strength in Fig. 1(c). For small interactions, 0 < g d < 0.002, the system is in a condensed SF and a single natural orbital is macroscopically occupied (λ 1 ≈ N ). With increasing interactions 0.002 < g d < 2, the system fragments and multiple natural orbitals have a macroscopic occupation. The occupations change gradually as the interactions increases and for g d 2 the crystal state is reached and m = N natural orbitals are (almost) equally populated [Fig. 1(b)]. The equal importance of the natural orbitals is a hallmark of the manybody features of the crystal state; it implies that the density matrices of higher order, ρ(p) (p > 1) are not representable as products of first-order density matrices ρ(1) . Thus the determination of quantities such as the phase diagram and the full distribution function of the atom number operator (see below) require a many-body model and cannot be achieved through the classical model we present in the SI [88], Sec. S2.

We now define the crystal state order-parameter ∆, (see Ref. [52])

∆ = k λ k N 2 , (2)
where λ k is the k th natural occupation. The maximal value, ∆ = 1, is obtained for the SF , while the CS is identified by the minimum value, ∆ = 1 N . Hence, for the N = 8 system at hand ∆ → 0.125 characterizes the CS transition. Fig. 3(a) shows the plot of ∆ as a function of g d and V clearly displaying the transition from the SM I state to the CS state.

For small interactions g d and a shallow lattice depths V , the bosons are fully condensed into a SF and ∆ ≈ 1. With increasing g d and/or V , the SM I forms; fragmentation and -consequently -a reduction in the value of ∆. Further increase of g d > 1 decreases ∆ gradually towards its minimum value ∆ = 0.125 for all values of V marking the onset of the KCS. When ∆ reaches its minimum, the maximally (eightfold) fragmented CS is reached; the orderings of the CS that we analyze in the following are the KCS or the DW CS. Importantly, by analyzing ∆ alone, the KCS → DW CS transition cannot be identified.

To identify the KCS → DW CS, we use the population imbalance of even and odd sites defined as

Θ = 1 N e,o n o -n e , (3)
where o and e label odd and even lattice sites, respectively, and n o and n e their respective population. The even-odd imbalance Θ is maximal for the density modulation corresponding to the DW CS state. Fig. 3(b) shows Θ as a function of g d and V to characterize the KCS → DW CS transition. The small imbalance for small values of g d results from the localization of the atoms in the central wells. For larger g d , Θ becomes small since the SM I possesses a uniform density.

As g d increases further, an increase to maximum evenodd imbalance indicates the KCS → DW CS transition. The KCS → DW CS transition shows a stronger dependence on V compared to SM I → KCS transition. A shallower lattice potential favors the KCS → DW CS transition at lower g d while deeper lattices require larger g d for the transition, Fig. 3(b).

We now discuss an experimental protocol to detect all the emergent states and thereby the state diagram presented above using standard imaging [82][83][84][START_REF] Ketterle | Inguscio[END_REF][START_REF] Smith | [END_REF]. These socalled single-shot measurements correspond to a projective measurement of the wavefunction. Ideally, the images contain an instantaneous snapshot of the position of all N particles distributed according to the N -particle probability distribution |Ψ| 2 . Here, we compute a set of single-shot simulations from our MCTDHB-groundstate wavefunctions [72][73][74]87] and evaluate the full distribution functions of the particle number, i.e., we quantify the probability P n (x) to detect n particles at positions x [Fig. 3(c-e)].

In the delocalized SM I several particles can be de-tected in the same site with a significant probability: P n (x) are nonzero for n ≤ 2. When the KCS is reached at g d ≈ 1, the bosons become completely separated and localized resulting in P n≥2 (x) ≈ 0. The transition from the SM I via the KCS to the DW CS with increasing interaction strength g d is characterized unequivocally through the analysis of the full distribution functions P 0 (x),P 1 (x) and P 2 (x). While P 0 and P 1 exhibit the distribution patterns of the KCS and the DW CS, P 2 → 0 signals the bosons complete isolation in crystal states.

The results of Fig. 3(c-e) show a good agreement with the ones of Fig. 1(a) and Fig. 3(a-b), demonstrating the KCS and DW CS transition for the same values of g d .

The simultaneous presence of density-wave order in P 0/1 and isolation, P n≥2 → 0, can thus experimentally identify the KCS and the DW CS (See SI [88], Sec. S4 for P n≥3 ≈ 0). An alternative experimental protocol using the variance of single-shot measurements to quantify the order parameter ∆ [52,73] intertwined with a binning of them to quantify the even-odd-imbalance Θ is described in the SI [88], Sec. S5.

We have thus proposed two viable experimental protocols to detect the state diagram of dipolar bosons in lattices with incommensurate filling, the crystal state transition, and the plethora of emergent crystal orders.

BC gratefully acknowledges the financial support from Department of Science and Technology, Government of India under DST Inspire Faculty fellowship. AUJL acknowledges financial support by the Austrian Science Foundation (FWF) under grant No. P-32033-N32, and the Wiener Wissenschafts-und TechnologieFonds (WWTF) project No MA16-066 ("SEQUEX). CL acknowledges funding by the FWF under grant No. M-2653. Computation time on the HPC2013 cluster of the IIT Kanpur and the HazelHen Cray computer at the HLRS Stuttgart is gratefully acknowledged. the potential energy, i.e., E K ≈ E I < E P in Fig. S2 for

g d ∈ [0.3, 2]
. For larger interaction energies the KCS order becomes energetically unfavorable and the DW CS state emerges. This DW CS state gradually melts as the interactions increase; when the interaction energy E I eventually dominates over the kinetic and potential energies, E I > E P > E K , a "pure" crystal state CS is formed.

S3. POSSIBLE CRYSTAL STATES

In this section, we investigate a complementary model of classical dipolar particles in a lattice in order to highlight that there are several possible ways that dipolar particles may arrange their positions in a lattice as a function of the strength of the dipolar interactions between them.

The energy of classical dipoles in a lattice is given as

E class (x 1 , ..., x N ) = N i=1 V ol (x i ) + i<j V int (x i -x j ). (1)
Here, x k is the position of the k-th particle, V ol (x) = V sin 2 (κx) is the lattice potential energy, and

V int = g d |xi-xj | 3
+α is the dipolar interaction between particles i and j. As in the main text, we use α = 0.05. Since we consider the ground state of classical particles (at rest), there's no kinetic energy contribution. Consequently, the kinetic crystal state (KCS) observed in the main text for quantum particles, is not present in our classical results below. However, we are able to demonstrate that the density-wave crystal state (DW CS) is a

A. Classical model verification

We now discuss the results for the positions of classical dipolar particles obtained from minimizing the energy in Eq. (1). To verify the classical model, we begin our investigation with the same configuration, N = 8 particles in S = 5 sites, as discussed in the main text for the quantum case. In Fig. S3, we show the positions x 1 , ..., x 8 that minimize E class as a function of the interaction for three values of the lattice depth, V = 5, 8 and V = 15.

The positions of classical dipoles in a lattice, due to the absence of kinetic do not feature a kinetic crystal state like the one observed for the quantum particles in the main text. However, the DW CS is obtained also in the classical model with the same ordering of particles ( [11], 1, [11], 1, [11]) similarly to the quantum case in the main text. We infer that we can use our classical model to investigate the possible DW CS orders -also for other configurations of the number of particles N and the number of lattice sites S. In Fig. S7(a-b), we show the probability P n (x) to find n particles at position x for the same parameters as in Fig. (3)(c-e) of the main text for P n>2 (x). Unlike that of the momentum space, the x-space probability becomes zero for P n>2 (x). There is a very small region where P 3 (x) = 0. P n≥4 (x) = 0 for all x.

S5. DETECTION OF THE CRYSTAL DENSITY WAVE FROM THE VARIANCE AND IMBALANCE IN SINGLE-SHOT IMAGES

Alternative to the analysis of the full distribution functions as shown in the main text, the emergent phases and the phase diagram for the one-dimensional lattices, incommensurately filled with dipolar bosonic atoms, can also be experimentally detected using the single-shot variance in combination with the single-shot expectation of the imbalance parameter Θ (Eq. (3) in the main text).

To characterize the transition from the SM I state to the CS, we use the variance V of simulated single-shot measurements in position space i.e. deviations of each single-shot measurement from the mean-value of many single-shot samples (see Refs. [1,2] for the mathematical definition of V).

The variance strongly depends on the degree of localization, i.e., on how much the positions of the bosons fluctuate from image to image [1]. For increasing localization, the variance decreases. Since the CS corresponds to a maximal localization of the system, V attains its minimum for the CS. In order to detect the transition to the CDW state, we evaluate the population imbalance parameter Θ, Eq. ( 3) of the main text, by binning the single-shot images: we count for each single-shot the number of atoms detected

FIG. 1 .

 1 FIG. 1. (a) One-body density ρ(x) as a function of dipolar interaction strength g d for a lattice depth V = 8. For small interactions g d 0.8, in the SM I state, the density exhibits a five-fold structure. As g d increases, the density develops a twofold splitting in the doubly-occupied central three wells displaying the onset of SM I → KCS transition. For even larger g d , the density distribution transitions to a pattern with alternating single and double occupations. This signifies the DW transition and formation of DW CS state. See Supplementary Information (SI) [88], Sec. S1 for the momentum space one-body density. (b) Natural occupations (plotted cumulatively) as a function of interaction strength g d . For small g d , the SF fraction results in the dominance of λ1. With increasing g d , the fragmentation increases as several λ k>1 become significant. Beyond the SM I → KCS transition, the system is maximally fragmented with M = 8 orbitals populated almost equally. Both x-axis in logscale.

FIG. 2 .

 2 FIG. 2. The two-body densities ρ2(x1, x2). (a) g d = 0.0005: The pure SF shows a nearly uniform density distribution. (b) g d = 0.1: SM I state shows a diagonal depletion because of the localization of bosons. (c) g d = 1.5: At the KCS, a correlation hole develops; the probability to find two bosons at the same place vanishes. (d) g d = 10.0: The DW CS shows a completely split but non-uniform density signaling density-wave order.

FIG. 3 .

 3 FIG. 3. Characterization and detection of the KCS and the DW CS. (a) Crystal order parameter ∆ as a function of interaction strength g d and lattice depth V . The maximum ∆ = 1 correspond to the completely condensed SF . The SM I is revealed by intermediate values of ∆. The minimum value of ∆ = 1 N = 0.125 corresponds to a CS thereby characterizing its formation. The SM I → KCS transition boundary is indicated by the continuous line (∆ = 0.2) while the KCS → DW CS obtained from Θ [Panel (b)] is indicated by the dashed line. (b) Imbalance parameter Θ as a function of g d and V . Both SM I and KCS correspond to low values of Θ. The transition to the maximum value of Θ indicates the KCS → DW CS transition. The continuous line displays the transition boundary (Θ = 3.0) while the dashed line shows the SM I → KCS transition obtained from ∆ [Panel (a)]. (c-e) Full distribution functions Pn(x) as a function of the interaction strength g d for a barrier height of V = 8 evaluated from 10000 single-shot images. The structure in the probability to find zero particles P0(x) and one particle P1(x) resembles that of the one-body density [compare Fig. 1(a)]. The vanishing of P2(x), the probability of finding two particles at x, clearly exhibits the transition to the CS. At large interactions, g d 6, P0(x) and P1(x) exhibit a density wave pattern. This pattern along with vanishing P2(x) reveals the DW CS state. See SI [88], Sec. S4 for the complementary full distribution functions in momentum space. All x-axis in logscale.

  FIG. S2. The kinetic, potential, and interaction energies, EK , EP , and EI , respectively, as a function of interaction strength g d for N = 8, V = 8, and S = 5. While the kinetic and potential energy EK + EP dominate for small to intermediate interactions, the interaction energy EI clearly dominates for large g d .

  FIG. S3. Verification of classical description. Plotted here, the classical position of the atoms x as a function of the interaction strength g d for different values of lattice depth V . The classical description displays different crystal arrangements: the DW CS state and the transition to the CS for strong interactions. Note the absence of the KCS state as well as the SM I state that is seen for the quantum case in the main text Fig. 1.

  FIG. S6. Full distribution functions Pn(k) in momentum space for the same parameters as in Fig. (3)(c-e) of the main text.

FIG

  FIG. S7. (a-b) Higher order distribution functions in real space Pn>2(x) for the same parameters as in Fig. (3)(c-e) of the main text. (c) Single shot variance as a function of position V(x) for V = 8.

  Fig. S7(c) displays the variance V as a function of position x for V = 8. A clear reduction of the variance V(x) is seen with interaction g d as the system transition from the SM I to the CS. The variance integrated over space V hence can be used to determine the CS transition experimentally. Fig. S8(a) displays the variance V as a function of g d and V . The SF state is maximally delocalized, thus having a maximum value of V. With increasing interactions, the SM I state is reached and the partial localization of the bosons at the lattice sites leads to a decrease of V. When the CS is reached at large values of g d , the bosons are completely localized. The vanishing overlap between the bosons forces V towards its minimum in the CS. The results of Fig. S8(a) are in good agreement with the ones of Fig. 3(a) of the main text, predicting the transition to the CS state for the similar values of g d and V . The variance of single-shot images can thus be used to experimentally determine accurately the transition to the crystal state.

Supplemental Material: Unveiling Emergent Crystal Orders of Incommensurate Dipolar Bosons in

One-Dimensional Lattices using Full Distribution Functions This Supplementary Information discusses the onebody momentum density in Sec. S1, the kinetic, potential, and interaction energies in Sec. S2, illustrates complementary possible crystal orders for different lattice sizes and particle numbers with a classical model in Sec. S3, as well as the higher-order distribution functions of the particle number operator in real space, P n>2 (x), and the full distribution functions of the particle number operator in momentum space in Sec. S4. An alternate detection protocol for the phase diagram [Fig. 3(a-b) of the main text] using the variance of single-shot images intertwined with a binning analysis is demonstrated in Sec. S5.

S1. ONE-BODY MOMENTUM DENSITY

Here, we analyze the one-particle momentum density for the same parameters as in Fig. 1 of the main text to assess the spatial coherence of the state, see Fig. S1.

For small to moderate interaction strengths (g d 0.8), the density exhibits a central momentum peak surrounded by a delocalized background distribution or small-amplitude peaks in a stripe-like topology. The amplitude and the width of the central peak reduce when the interaction strength g d increases. The presence of a localized central peak indicates the presence of the delocalized SF fraction that coexists with the localized M I state which in turn is responsible for the delocalized background distribution. When the transition to the crystal state occurs at g d ≈ 1, the momentum distribution completely loses the stripe-like topology forming a uniformly delocalized distribution: the crystal state is character- * bchat@iitk.ac.in † Schmiedmayer@atomchip.org ‡ camille.leveque@tuwien.ac.at § axel.lode@univie.ac.at ized by strong spatial confinement and decoherence, i.e., a vanishing superfluid fraction.

S2. ENERGY

The emergence of distinct quantum phases in dipolar interacting ultracold bosons is a result of the interplay between the interaction, kinetic, and potential energies. In this section we discuss this interplay of energies for the same system as shown in the main text, i.e., N = 8, V = 8, S = 5, as a function of the interaction strength, see Fig. S2.

The kinetic and potential energy clearly dominate over the interaction energy at small to intermediate interaction strength g d 0.3. This dominance is responsible for the formation of the SF and SM I states. The KCS state forms in the region, where the interaction energy and the potential energy are comparable, but smaller than

B. Density-wave crystal states for different filling factors

We now discuss the density-wave crystal states that we obtain from our classical model, Eq. ( 1), for the same number of lattice sites as in the main text, S = 5, but for different incommensurate filling factors, see Fig. S4.

The DW CS states with orders [11], [1], [11], [1], [11] and , and (c), respectively. Note that, for a larger particle number, the transition from a DW CS to a CS occurs for lower values of interaction g d . This is because interaction energy per site grows faster compared to the potential and kinetic energies as a function of particle number at a fixed value of g d .

C. Density-wave crystal states for different lattice sizes

We now discuss the density-wave crystal states that we obtain from our classical model, Eq. ( 1), for larger lattice sizes than in the main text (S = 7, 9), but for a similar incommensurate filling factor, see Fig. S5.

It is clearly seen that the DW CS with the configuration ..., [11], 1, [11], 1, [11], ... prevails also for the case of lattices with a larger number of sites.

In summary, we find that there is a rich variety of density-wave crystal states of dipolar particles, already in the classical model that we investigated here. This underlines the importance of the experimental protocols that we put forward to detect these states in quantum systems.

S4. FULL DISTRIBUTION FUNCTIONS IN MOMENTUM SPACE Pn(k) AND IN REAL SPACE FOR LARGER Pn>2

A.

Full distribution functions in momentum space

In Fig. S6, we show the probability P n (k) to find n particles with momentum k for the same parameters as in Fig. (3)(c-e) of the main text. P 0 (k) and P 1 (k) shows identical distribution with the momentum density [Fig. 1(b)] of the main text but with opposite signs.

Increasing n shows a distinct narrowing of the distribution. For higher-order distributions P n>3 (k), the maximum probability is centered at k = 0 but falls off sharply and the background is strongly reduced. We remark here, that the total number of events in the plots for P 5 (k) and, especially, for P n≥6 (k) is extremely small; thus the probability to find 5 or more particles at the same momentum is extreme small. For P 5 (k) the maximum count is O(10) and for P 6 (k) it is O(1) out of N × N shots = 8 × 10000 = 80000.