# Correlations of strongly interacting one-dimensional ultracold dipolar few-boson systems in optical lattices

Budhaditya Chatterjee, Marios Tsatsos, Axel Lode

► To cite this version:

Budhaditya Chatterjee, Marios Tsatsos, Axel Lode. Correlations of strongly interacting one-dimensional ultracold dipolar few-boson systems in optical lattices. New Journal of Physics, 2019, 21, pp.033030. 10.1088/1367-2630/aafa93. hal-02370334

HAL Id: hal-02370334

https://hal.science/hal-02370334

Submitted on 19 Nov 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Correlations of strongly interacting one-dimensional ultracold dipolar few-boson systems in optical lattices

To cite this article: Budhaditya Chatterjee et al 2019 *New J. Phys.* 21 033030

View the [article online](#) for updates and enhancements.

## Recent citations

- Management of the correlations of UltracoldBosons in triple wells
  Sunayana Dutta et al

- Many-body quantum dynamics of an asymmetric bosonic Josephson junction
  Sudip Kumar Haldar and Ofir E Alon
New Journal of Physics

The open access journal at the forefront of physics

Deutsche Physikalische Gesellschaft DPG

IOP Institute of Physics

Published in partnership
with: Deutsche Physikalische
Gesellschaft and the Institute
of Physics

<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED
15 June 2018

REVISED
31 October 2018

ACCEPTED FOR PUBLICATION
21 December 2018

PUBLISHED
28 March 2019

Original content from this
work may be used under
the terms of the *Creative
Commons Attribution 3.0
licence*.

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

# Correlations of strongly interacting one-dimensional ultracold dipolar few-boson systems in optical lattices

Budhaditya Chatterjee¹,⁵, Marios C Tsatsos² and Axel U J Lode³,⁴

¹ Department of Physics, Indian Institute of Technology-Kanpur, Kanpur 208016, India

² São Carlos Institute of Physics, University of São Paulo, PO Box 369, 13560-970, São Carlos, SP, Brazil

³ Wolfgang Pauli Institute c/o Faculty of Mathematics, University of Vienna, Oskar-Morgenstern Platz 1, A-1090 Vienna, Austria

⁴ Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee 2, A-1020 Vienna, Austria

⁵ Author to whom any correspondence should be addressed.

E-mail: bchat@iitk.ac.in, marios@ifsc.usp.br and axel.lode@univie.ac.at

**Keywords:** dipolar boson, quantum phases, coherence properties, correlation functions, optical lattice

Abstract

Strongly interacting finite ensembles of dipolar bosons in commensurately filled one-dimensional optical lattices exhibit diverse quantum phases that are rich in physics. As the strength of the long-range boson-boson interaction increases, the system transitions across different phases: from a superfluid, through a Mott-insulator and a Tonks-Girardeau gas to a crystal state. The signature of these phases and their transitions can be unequivocally identified by an experimentally detectable order parameter, recently described in *Phys. Rev. A* **98** 235301 (2018 [33]). Herein, we calculate the momentum distributions and the normalized Glauber correlation functions of dipolar bosons in a one-dimensional optical lattice in order to characterize all their phases. To understand the behavior of the correlations across the phase transitions, we first investigate the eigenfunctions and eigenvalues of the one-body reduced density matrix as a function of the dipolar interaction strength. We then analyze the real- and momentum-space Glauber correlation functions, thereby gaining a spatially and momentum-resolved insight into the coherence properties of these quantum phases. We find an intriguing structure of non-local correlations that, independently of other observables, reveal the phase transitions of the system. In particular, spatial localization and momentum delocalization accompany the formation of correlated islands in the density as interactions become stronger. Our study showcases that precise control of intersite correlations is possible through the manipulation of the depth of the lattice, while intrasite correlations can be influenced by changing the dipolar interaction strength.

**1. Introduction**

Ultracold atoms with dipole-dipole interactions have become a popular tool to simulate and understand the physics of long-range interacting systems [1, 2]. The experimental realization of quantum gases having a dipole-dipole interaction has been achieved with atoms having permanent magnetic dipole moments, such as chromium [3, 4], dysprosium [5] and erbium [6] as well as with polar molecules, for instance, potassium-rubidium [7] and cesium-rubidium [8]. Owing to the long-range and anisotropic nature of dipole-dipole interactions, novel quantum effects not present in atoms with contact interactions emerge. Prominent examples include the elongation of the condensate along the orientation of the dipole moments [9–11] and the exciting phenomenon of geometrical stabilization [2] of a dipolar Bose-Einstein condensate in traps of certain shapes, like extremely oblate ones [9, 10, 12–14]. More recently, the fascinating formation of quantum droplets in dipolar condensates has been observed and then theoretically studied [15–18].

In ultracold systems, the dimensionality is an experimental control-parameter and crucially important:
lower-dimensional systems often produce a variety of effects not seen in three spatial dimensions. The
occurrence of *p*-wave superfluidity in two-dimensional Fermi gases [19, 20] provides such an example. Dipolar
atoms in quasi-one-dimensional traps are more amenable experimentally since the collisional instabilities arising from the head-to-tail alignment in two and three spatial dimensions are prevented [14, 21]. Unidimensional dipolar atoms have been predicted to exhibit Luttinger liquid-like behavior [22–25] as well as anisotropic effects in curved and ring geometries [26–28]. Moreover, for very strong dipolar interactions a remarkable crystallization effect takes place where the dipolar atoms themselves form a crystal lattice structure irrespective of the geometry of their external confinements [22, 26, 29–33].

Optical lattices often serve as a controllable toolbox to understand and simulate a large variety of condensed matter systems. For dipolar atoms, the additional existence of the long-range anisotropic interactions leads to a plethora of interesting quantum phases arising from the interplay of the kinetic energy, the short and long-range interactions, each dominating different energy scales [2]. Density waves [34, 35], Haldane insulators [35, 36], checkerboard patterns [34, 37] and Mott solids [38] are some prominent examples of these phases.

In our study, we consider a system where four different phases are amalgamated: superfluid (SF), Mott insulator [39], fermionized Tonks–Girardeau gas [40–43] and a crystal-like state [26, 31, 43] can each emerge in a finely tuned system of few dipolar bosons in a multi-well trap. Superfluidity appears due to the bosonic nature of the particles combined with their weak interactions. A finite well depth together with somewhat stronger interactions breaks superfluidity and leads the system to behave as a so-called Mott-insulator (MI) [44, 45]. For moderately strong interactions, the system mimics the boson-to-fermion mapping, that is known to apply exactly at the infinite-strength limit of contact interactions as well as finite-size hard-rod interactions [46]. In this Tonks–Girardeau limit, the particles isolate themselves from their neighbors in order to avoid infinities in the interaction energy. Thus the bosonic density approaches that of its non-interacting fermionic counterpart. Nevertheless, the momentum density of the bosonic system is still distinct [46]. We remark here that an exact Bose–Fermi map for long-range dipolar interaction also exists, exploiting the divergence of the dipolar interaction at zero separation [30, 31]. This dipolar Bose-Fermi map implies that bosons with strong dipolar interactions and fermions with strong dipolar interactions have identical one-body densities and identical one-body momentum densities. Last, for even higher interaction strengths, the long-range tail of the interaction dominates and leads to the formation of the so-called crystal phase [22, 26, 29–33].

Herein we follow the same strategy as [47–49] and theoretically investigate, the physics of a larger many-body system by studying in detail its few-body building blocks. Although few-body systems with dipole–dipole interactions have not been studied experimentally as of yet, the experimental realization of tunable few-fermion systems [47–49] motivates our theoretical study of a small ensemble of dipolar bosons. Moreover, our results for the few-atom system can be generalized to larger systems (see appendix A). We investigate a system of dipolar bosons in an optical lattice by studying the triple-well potential. Theoretically, dipolar atoms in triple-well traps have been explored using mean-field methods [50], the extended-Bose–Hubbard model (BHM) [51–56] and also the multi-configurational time-dependent Hartree (MCTDH) method [32]. Notably, the mean-field methods and BHM are unable to address very strong dipolar interactions. It is thus necessary to employ a general many-body approach for the cases where the strong dipolar interactions dictate the physics of the system.

The MCTDH for bosons (MCTDHB) [57, 105] is such a general many-body method capable of addressing strong interaction regimes [58–62] and its implementation in the MCTDH-X software [63–65] has been employed in [33] to establish an order parameter and an experimental method to classify and detect all the quantum phases of dipolar atoms in optical lattices via single-shot images.

In this paper, we explore the transitions across the above-mentioned four phases by varying the interaction strength and then analyze the normalized Glauber correlation functions to unveil the coherence properties of each phase. We show how the system can be brought to any desired phase and how the transitions are reflected in the correlation landscape in coordinate and momentum space. We observe distinct structural changes in the correlations that accompany the phase transitions and clearly characterize the phases and their transitions. We stress here that the focus of this work is different from that of [33]. While in [33], an order parameter and its experimental detection via single-shot images was the prime motivation, here we take the emergence of these phases as the starting point and focus on their correlation and coherence properties in coordinate and momentum space.

We remark that the few-body finite-size system studied here cannot exhibit true phase transitions in the macroscopic sense. What we obtain are ground states that are the finite-size precursors to the macroscopic thermodynamic phases [66]. For the sake of simplicity and as an analogy to the thermodynamic systems, we still use the term ‘phase’ to refer to them.

This paper is structured as follows: section 2 introduces the Hamiltonian and quantities of our interest. Section 3 discusses the eigenfunctions and eigenvalues of the reduced one-body density matrix for the emergent phases. Section 4 presents an analysis of the normalized Glauber correlation functions. Section 5 discusses the effect of lattice depth and section 6 provides an outlook concluding our work.
## 2. Model

Our system consists of N polarized, dipolar bosons of mass M in a one-dimensional optical lattice and is governed by the Hamiltonian

$$H = -\sum_{i=1}^{N} \frac{\hbar^2}{2M} \partial_{x_i}^2 + \sum_{i=1}^{N} V_{\text{ol}}(x_i) + \sum_{i<j} V_{\text{int}}(x_i - x_j). \quad (1)$$

The one-body potential $V_{\text{ol}}$ represents a quasi-one-dimensional optical lattice potential, modeled as $V_{\text{ol}} = V \sin^2(\kappa x)$. Here, V is the depth of the lattice and $\kappa$ its wave number. In order to confine the bosons to a desired number of sites, we impose a hard wall boundary condition at $x = \pm S\pi/2\kappa$ where S is the number of lattice sites (for odd S). A strong transverse confinement of characteristic length $a_\perp$ ensures that the system is quasi-one-dimensional by preventing excitations in the transverse direction.

We consider a pure dipole-dipole interaction where the interaction potential can be written as $V_{\text{int}}(x_i - x_j) = \frac{g_d}{|x_i - x_j|^3 + \alpha}$. For large separations $|x_i - x_j| \gg a_\perp$, the interaction potential varies as $V_{\text{int}}(x_i - x_j) \sim 1/|x_i - x_j|^3$. For small separations $|x_i - x_j| \lesssim a_\perp$, the transverse confinement induces a short-scale interaction cutoff $\alpha \approx a_\perp^3$, thus regularizing the divergence at $x_i = x_j$ [31, 67, 68]. Here $g_d$ is the dipolar interaction strength, given as $g_d = d_m^2/4\pi\epsilon_0$ for electric dipoles and as $g_d = d_m^2\mu_0/4\pi$ for magnetic dipoles, $d_m$ being the dipole moment, $\epsilon_0$ the vacuum permittivity and $\mu_0$ the vacuum permeability. In order to obtain universal dimensionless quantities we rescale equation (1) in terms of the recoil energy $E_R = \hbar^2\kappa^2/2M$ effectively setting $\hbar = M = \kappa = 1$. All quantities become dimensionless with the length given in units of $\kappa^{-1}$.

We remark that, generally, the dipole-dipole interaction potential in 1D also includes a contact (Dirac delta) term owing to the transverse confinement [31, 67, 68]. However, in order to examine the effect of the finite-range interaction alone, we set the contact interaction term equal to zero. In an experiment, this could be achieved via management of Feshbach resonances.

We obtain the ground states of the dimensionless Hamiltonian using the imaginary time-propagation of the MCTDHB equations of motion [57]. For the system studied, seven orbitals are sufficient to obtain converged solutions for the superfluid, Mott-insulating, and crystal phases. Our motivation is to understand the many-body physics from finite size few-body elemental building blocks and so it is necessary to consider a minimal system that displays all the fundamental physics. We thus consider a triple well setup ($S = 3$) with $N = 6$ since it represents the minimal setup to understand all the relevant physics for the dipolar boson in an optical lattice. Similarly, we choose the barrier height $V = 8$ and short-range cutoff $\alpha = 0.05$ such that the emergent physics are exhibited distinctly. A setup with these values is within experimental feasibility. Our results are fundamental, and they can be readily generalized for larger systems, as we demonstrate in appendices A, B, and C.

## 3. Natural orbitals, natural occupations and quantum phases

While, in principle, the wavefunction contains all information about a quantum system, the aspects of correlations are better assessed through the reduced density matrix (RDM). Considering an N-particle state $|\Psi\rangle$, the pth order RDM is obtained by tracing $N - p$ degrees of freedom [69, 70] as:

$$\rho^{(p)} = \text{tr}_{p+1 \dots N} [ |\Psi\rangle \langle\Psi| ]. \quad (2)$$

For $p = 1$, we obtain the 1st order or one-body RDM

$$\rho^{(1)}(x, x') = \langle \hat{\Psi}^{\dagger}(x) \hat{\Psi}(x') \rangle = \int \Psi(x_1, \dots, x_N) \times \Psi^*(x'1, \dots, x_N) dx_2 \dots dx_N, \quad (3)$$

where $\Psi \equiv \Psi(x_1, \dots, x_N)$ is the many-body wavefunction and $\hat{\Psi}^{(\dagger)}$ the field annihilations (creation) operator in 1D position space.

In the following, we analyze the one-body RDM expanded in its eigenfunctions $\varphi_i$ as:

$$\rho^{(1)}(x, x') = \sum_{i} \lambda_{i} \varphi_{i}^{*}(x) \varphi_{i}(x'). \quad (4)$$

The eigenfunctions $\varphi_i$ are called natural orbitals and the corresponding eigenvalues $\lambda_i$ are called natural occupations. Each $\lambda_i$ represents the population of the ith orbital. The spectral decomposition of the one-body RDM is particularly useful since it serves to define Bose-Einstein condensation (BEC) in an interacting many-body system: if the largest natural occupation is of the order of the number of particles N, the system is said to be condensed [71]. If there are more than one natural orbitals with populations of order N then the system is said to be fragmented [72-74, 106-111].
<figure><img src="image_4.png" /></figure>

The natural orbitals and their populations are very important in characterizing the phases occurring for dipolar bosons in an optical lattice. The natural occupations directly relate to the emerging phases of strongly interacting dipolar bosons [33, 61, 75]. Figure 1 displays the evolution of the natural occupation as a function of interaction strength $g_d$.

In a conventional BEC with contact interactions, the interplay between contact interactions, potential and kinetic energy determines the quantum properties [44, 45]. For dipolar bosons, the additional long-range interactions make it a four-way competition between the kinetic and potential energy, the short-range and the long-range interactions that lead to the existence of various phases at different energy scales [33, 76]. When the interaction is very small ($g_d \approx 0$) the kinetic energy dominates, thus leading to a SF phase. Here, the system forms a BEC since only the first natural orbital is macroscopically populated. Hence, in this region $\lambda_1 \approx N$. As $g_d$ increases the interaction energy starts to dominate over the kinetic energy and the system arrives at the MI phase [39, 45, 61, 77–79]. In the MI phase, the bosons are localized inside each well and occupy S (i.e. a number equal to
the number of sites) orbitals equally, hence $\lambda_i = N/S$. At the same time tunneling between the lattice sites is strongly reduced [33]. As $g_d$ increases further, the long-range effect of the interaction now affects the many-body state and the crystal phase is reached [26, 31, 43]. The already localized bosons in each well maximize their spatial separation due to the strong repulsion. The one-body RDM shows maximal fragmentation with $N$ orbitals equally populated, i.e. $\lambda_i = 1$, $i = 1, ..., N$. The formation of the crystal phase is driven only by the long-range interaction potential and is independent of the lattice potential $V_{\text{ol}}$.

Note that the TG gas acts as an intermediate step in the MI → CS transition. Unlike the SF, MI and CS states, a universal behavior of the natural occupation does not exist for the TG gas. Here, the highest natural occupation $\lambda_1$ scales as $N^{-1/2}$ (N being the particle number) [81]. Thus, the formation of the TG gas cannot be determined by the natural occupations alone. Nevertheless, the TG gas can be distinguished by its density and correlation signatures. Interestingly, for a vanishing cutoff $\alpha \to 0$, one cannot obtain an interacting SF state or an MI state. When $\alpha \to 0$, the interaction diverges at $x_i = x_j$, thus any non-zero $g_d$ will result in the system being either a TG gas or a super-Tonks Gas depending on the interaction sign and density [30, 31, 82]. Since we have a non-zero cutoff $\alpha$, the divergence at $x_i = x_j$ is avoided and we can get both SF and MI phases by suitably changing $g_d$. We remark here that Bloch oscillations in the absence of a lattice have recently been observed in a Tonks gas [80].

We remark that the BHM cannot describe the crystal phase. The BHM employs a tight-binding approximation to construct a site-localized Wannier basis. The BHM cannot describe structures smaller than the lattice spacing since the localized on-site Wannier basis are single-particle functions that do not show the nodal structures that result from strong dipolar interactions. The single-band BHM cannot describe the physics when the energy becomes comparable to the higher bands [83]. The crystal phase occurs for very strong dipolar interactions where the contribution of the higher bands cannot be neglected and the many-body wavefunctions possess structures smaller than the lattice spacing. Thus the single-band BHM cannot describe the crystal phase and its transition.

Moreover, for strong dipolar interactions, the interaction energy dominates the potential energy; the crystal formed by the bosons can therefore have a 'lattice constant' that depends on the interaction strength. The physics of such a crystal cannot be cast into a model where the basis is localized in space and the corresponding lattice constant is fixed irrespective of the interaction strength, as done in the single-band and multi-band BHM. To represent the many-body state using a multi-band BHM for a system where the lattice sites are not aligned with the actual positions of the particles, a very large number of bands becomes necessary quickly. Thus, even the multi-band BHM—albeit a formally exact description when an infinite number of bands is allowed—is not the method of choice for the crystal state at large dipole-dipole interaction strengths.

The following schematic picture qualitatively illustrates the basic mechanism of fragmentation and fermionization for six particles in a triple well. At the non-interacting limit, the ground state is in a single superposition $\phi_{GP}$ of the three modes ($\phi_l, \phi_c, \phi_r$) that are localized in the left, center and right wells, respectively. The existence of hard walls at the border of our potential breaks the translational invariance of the Hamiltonian and so the central well is slightly denser compared to the other two. Neglecting this asymmetry, the particle configurations (i.e. distributions of bosons across the wells) are degenerate energetically. The bosonic superposition $\phi_{GP}$ is coherent and the ground state is condensed and SF. In other words, a single orbital $\phi_{GP} = \frac{1}{\sqrt{3}}(\phi_l + \phi_c + \phi_r)$ is enough to describe the ground state of the many-body system, as the product state $\prod_{i=1}^{N} \phi_{GP}(x_i) = |N, 0, 0, ..., 0\rangle$. Now consider a small but finite interaction parameter $g_d$. If the lattice sites are firmly separated from each other via sufficiently high barriers then the orbitals $\phi_l, \phi_c, \phi_r$ are not overlapping. In that case, the threefold fragmented state $\Psi = |\frac{N}{3}, \frac{N}{3}, \frac{N}{3}, 0...0\rangle$ is energetically favorable for any finite interaction strength against $\Psi_{GP} = |N, 0, 0, ..., 0\rangle$ since the fragmented configuration minimizes the interaction energy (see also [84]). Hence, for large barrier heights $V$, a tiny interaction among the particles disrupts condensation in favor of fragmentation. The threefold shape of the trap determines that the number of significantly occupied orbitals is three. For a total number of $N = 6$ particles, two particles occupy each lattice site of the trap. For a further growth of the interaction, each of the two particles in the same site starts to repel its neighbor considerably; the orbitals now tend to avoid overlap in order to minimize the (otherwise increasing) interaction energy. This is the onset of fermionization, i.e. a density profile where bosons mimick a non-interacting fermionic density. For higher $g_d$, where the long-range tail of the interaction potential begins to dominate, we enter the crystal phase where the orbitals completely dissociate and each boson occupies a separate orbital, giving rise to a six-fold fragmented state. Here, the particle distribution and orbital shape are determined by the long-range interaction solely and not the lattice potential.

The above picture with single-configurational states is certainly oversimplified since the interactions change the shape of the natural orbitals along with their occupations (figure 1(b)). For finite interactions, there is a substantial intermixing of the non-interacting orbitals, which is parity preserving and thus a numerically exact
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a)–(d): The 1st order spatial correlation function $|g^{(1)}(x_1, x'_1)|$ and (e)–(h): the 2nd order correlation function $|g^{(2)}(x_1, x_2)|$, shown for all emergent phases. At $g_d = 0$, the bosons are coherent, a fact reflected in the uniform distributions: $|g^{(1)}| = 1$ and $|g^{(2)}| = 1$ everywhere. For the Mott insulator ($g_d = 0.1$), the bosons are completely localized in each well and thus the long-range off-diagonal coherence is destroyed: first-order coherence is restricted to individual lattice sites while second-order coherence is characterized by particle bunching on the off-diagonal, i.e. $|g^{(2)}(x_1, x_2)| > 1$ where $x_1$ and $x_2$ are positions in different sites and anti-bunching within individual lattice sites, i.e. $|g^{(2)}(x_1, x_2)| < 1$ for $x_1 \sim x_2$. For larger values of $g_d$ the intrawell first-order coherence decreases further. Also, the anti-bunching within sites is augmented ($|g^{(2)}(x_1, x_2)| \to 0$ for $x_1 \sim x_2$) while the bunching between distinct sites disappears as the interaction strength increases ($|g^{(2)}(x_1, x_2)| \to 1$ for positions $x_1$ and $x_2$ in distinct sites). All units shown are dimensionless.</figcaption></figure>

solution must be sought. However, the above line of thought illustrates the idea of the emergence of the different phases and their relation to fragmentation.

In figure 1(b) we plot the evolution of the natural orbitals $\phi_i(x)$ scaled by the natural occupations $\lambda_i$ with the interaction strength $g_d$. The orbital $\phi_1$ exhibits the mean-field condensate physics at small interaction strengths. For weak interactions, it shows the larger center-well population expected in the SF phase. However, it is unable to predict the equalizing of populations across the three wells in the MI and crystal phases. While $\phi_1$ alone captures the mean-field perspective, all three lowest orbitals $\phi_{1,2,3}$ are necessary to account for the SF-MI transition and, thus, to understand the single-band Bose-Hubbard picture. In other words, the system fragments with equal natural occupations $\lambda_1, \lambda_2, \lambda_3$ across the three natural orbitals whose density maxima are located at the minima of each well. For strong $g_d$, $\phi_2$ and $\phi_3$ show a splitting of the outer maxima. The higher orbitals $\phi_{4,5,6}$ contribute significantly only at the strongly interacting crystal regime. The primary feature of the latter is the enhancement of the intrawell (i.e. referring to the same well) splitting at strong interactions. Most importantly, the crystal phase is realized by the equal contribution of N natural orbitals [33].

# 4. Correlation functions and quantum phases

We now discuss the many-body coherence and correlation properties of the system in its distinct phases. To this end, we assess the coordinate and momentum space correlations in the state $|\Psi\rangle$, by studying the 1st and 2nd order Glauber correlation functions [70].

## 4.1. Spatial 1st order correlations

The 1st order Glauber correlation function

$$g^{(1)}(x'_1, x_1) = \frac{\rho^{(1)}(x_1, x'_1)}{\sqrt{\rho(x_1)\rho(x'_1)}}, \qquad (5)$$

is constructed by normalizing the one-body RDM (equation (3)) to a product of the respective one-body densities. It measures the proximity of the many-body state to a mean-field state in which $\rho^{(1)}(x_1, x'_1) = \rho^{(1)}(x_1, x'_1) = \rho(x_1)$ (i.e. uniform anti-diagonals) [61, 64, 65, 69, 85–87]. $g^{(1)}$ is generally a complex quantity (for ground states it is typically real) and is associated with phase coherence, which can be accessed
through interference experiments [88–91]. $g^{(1)}$ is bounded within [0, 1], with $|g^{(1)}| = 1$ implying perfect coherence and $|g^{(1)}| = 0$ the complete absence of it. By construction, the diagonal of $|g^{(1)}|$ is always equal to unity, so the coherence investigation needs to be done on the off-diagonal: 1st order coherence between any two points $x_1, x'_1$ denotes off-diagonal long-range order (ODLRO) [104].

Figure 2 displays $|g^{(1)}(x_1, x'_1)|$ for various interaction strengths that correspond to the emerging distinct phases. For $g_d = 0$ the bosons are delocalized over the whole lattice and are fully coherent. The many-body state is exactly described by a mean-field product of a single orbital $\phi_{GP}$ (see also discussion in section 3). The complete coherence is reflected in the uniform distribution of $|g^{(1)}(x_1, x'_1)| = 1$ throughout the lattice. As the interaction is switched on, the ODLRO is disrupted, thus reflecting the particle localization at each well. In an MI state at $g_d \approx 0.1$, the bosons are completely localized and, as a consequence, coherence is strongly reduced; $|g^{(1)}|$ now shows a strong off-diagonal reduction. However, the localization is still partial, since inside each well, the two residing atoms are delocalized and exhibit relative coherence. This intrasite delocalization and coherence is seen in the block-diagonal form of $|g^{(1)}|$ in figure 2(b). Here, each coherent block, occupying a space of roughly $\frac{\pi}{2} \times \frac{\pi}{2}$ square units, pertains to each lattice site.

As the interaction strength $g_d$ increases, its short-range portion dominates$^6$, leading to a fermionization-like behavior of the bosons [40–43] and a reduction of the intrawell coherence. At $g_d = 1.0$, which falls in this Tonks–Girardeau-gas-alike regime, one can already see the off-diagonal contribution vanishing (figure 2(c)).

For a further rise of the interaction strength $g_d$, the long-range tail of the interaction begins to dictate the physics. The influence of the strong long-range tail of the interaction potential forces the localized bosons maximally apart, thus avoiding each other. As a consequence, the coherence decreases even further.

Last, the system enters the crystal phase at $g_d = 15$, where complete fragmentation characterizes the state. The bosons are entirely localized even within wells, and the intrawell coherence is zero. The coherence blocks of $|g^{(1)}|$ are now centered at each boson location (figure 2(d)). At the crystal phase, the bosons maximize their interparticle distance and do no longer necessarily obey the external lattice spacing. This independence of the lattice spacing is seen, for instance, in the reduced size of the outermost blocks on the diagonal of $g^{(1)}$ ($x_1 \approx x'_1$) where $g^{(1)} \approx 1$ in figure 2(d) as compared to figure 2(c). The absence of off-diagonal coherence reflects the complete coordinate-space localization. Exceptionally, coherence marginally remains at a small region around the diagonal.

## 4.2. Spatial 2nd order correlations

The second-order RDM $\rho^{(2)}(x_1, x_2, x'_1, x'_2)$ quantifies the correlations between two particles, with its diagonal kernel $\rho^{(2)}(x_1, x_2) = \rho^{(2)}(x_1, x_2, x'_1 = x_1, x'_2 = x_2)$ representing the probability of the simultaneous detection of a particle at $x_1$ and another particle at $x_2$.

By normalizing $\rho^{(2)}(x_1, x_2, x'_1, x'_2)$ in terms of products of the respective one-body densities, we obtain the 2nd order Glauber correlation function

$$g^{(2)}(x_1, x_2, x'_1, x'_2) = \frac{\rho^{(2)}(x_1, x_2, x'_1, x'_2)}{\sqrt{\rho(x_1)\rho(x'_1)\rho(x_2)\rho(x'_2)}}, \quad (6)$$

which quantifies the 2nd order coherence in the system. Hereafter, we use the diagonal $|g^{(2)}| \equiv |g^{(2)}(x_1, x_2)| = |g^{(2)}(x_1, x_2, x'_1 = x_1, x'_2 = x_2)|$.

We would like to stress the fact that $\rho^{(2)}(x_1, x_2)$ and $g^{(2)}(x_1, x_2)$ provide distinct information about the physics of the state. $\rho^{(2)}(x_1, x_2)$ gives the spatial two-body probability distribution, but provides no information about the coherence. In contrast, $g^{(2)}(x_1, x_2)$, the second-order normalized Glauber correlation function, contains essential information about the coherence properties of the system that cannot be observed solely from $\rho^{(2)}(x_1, x_2)$; second-order coherence [$g^{(2)}(x_1, x_2) = 1$] implies that the two-body density is just a product of one-body densities, i.e. $\rho^{(2)}(x_1, x_2) = \rho^{(1)}(x_1, x_1) \rho^{(1)}(x_2, x_2)$, while second-order incoherence [$g^{(2)}(x_1, x_2) \neq 1$] implies that $\rho^{(2)}(x_1, x_2)$ contains information that cannot be extracted from $\rho^{(1)}$ alone. The discussion of the two quantities, $\rho^{(2)}(x_1, x_2)$ and $g^{(2)}(x_1, x_2)$, is thus independent. An analogous argument holds for the momentum space representation of the normalized Glauber correlation functions that we discuss in the following section.

For weak interactions ($g_d \approx 0$) the state is at the SF phase yielding $|g^{(2)}| = 1$, a fact that demonstrates 2nd order coherence and the absence of any correlation in the measurement of the positions of any pair of bosons, figure 2(e). In general, a SF phase shows Nth order coherence [69] and any p-particle detection probabilities are not correlated. As the interaction increases, the diagonal of $|g^{(2)}|$ displays regions with anticorrelations $|g^{(2)}| < 1$ (figure 2(f)) thus demonstrating loss of second-order coherence and localization, or else anti-bunching. In analogy to photon count statistics (see [93] and references therein), we use here the term (anti-)bunching to

<footnote>$^6$The ‘short-range part of the dipolar interaction’ refers to the part of the interaction potential within a single well; quantified as, $V_S \equiv \int_0^{\frac{\pi}{2}} \frac{g_d}{x^3+\alpha}dx$. Similarly, the ‘long-range part of the interaction’ refers to the part of the interaction potential beyond the local well; quantified as, $V_L \equiv \int_{\frac{\pi}{2}}^{\infty} \frac{g_d}{x^3+\alpha}dx$.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a)–(d): 1st order momentum correlations and (e)–(h): 2nd order momentum correlations for the distinct emergent phases of dipolar bosons in lattices. In the Mott-insulating state $g^{(1)}(k'_1, k_1)$ shows a clean diagonal pattern of stripes with alternating weak $|g^{(1)}| \approx 1$ and strong $|g^{(1)}| \approx 0$ correlations. As the interaction strength is increased to $g_d = 1.0$, the diagonal correlation hole (black regions with $|g^{(1)}| \approx 0$) spreads. A striped pattern similar to the Mott-insulating phase, but with a different spacing is seen in the crystal phase at $g_d = 15$. The 2nd order momentum correlations $g^{(2)}(k_1, k_2)$ show similar patterns, but with their value lying in the interval [0, 2] thus signifying (anti-)correlation, i.e. (anti-)bunching for $g^{(2)} > 1$ ($g^{(2)} < 1$).</figcaption></figure>

denote (decreased) increased value of $g^{(2)}(x_1, x_2)$, as compared to unity. In the context of quantum optics, (anti) correlated emission events determine the detection probabilities and their distributions. Similarly, two bosons will bunch together and, hence, localize if $|g^{(2)}| > 1$. At the MI state, an anti-bunching block covers each lattice site, marking lack of second-order coherence between the localized bosons in each site. For stronger interactions, the localized bosons within a site form intrawell correlation holes (see figure 2(g)). These intrawell structures are the key signature of higher-band effects that manifest both in the intermediate fermionization effect at $g_d = 1$ (figure 2(g)) and the emergence of the crystal phase at $g_d = 15$ (figure 2(h)). While the diagonal pattern of $|g^{(2)}|$ shows no significant differences between the $g_d = 1.0$ and the $g_d = 15$ cases, the off-diagonal shows a significant reduction at $g_d = 15$. Furthermore, the outermost antibunching blocks on the diagonal of $g^{(2)}$ are reduced in size and all blocks no longer aligned with the minima of the external lattice potential. Thus, the second-order correlations in figures 2(g), (h) display the emergence of extreme localization and dominating long-ranged interactions characteristic for the crystal phase.

### 4.3. Momentum correlations

We complement our investigations of the spatial correlations by analyzing the correlations $\tilde{g}^{(1,2)}(k, k')$ of the state $\Psi$ in momentum space. $\tilde{g}^{(p)}(k_1, \dots)$ is found from the reduced densities of the Fourier transformed many-body state $\tilde{\Psi}(k_1, \dots)$. We remark that the diagonals of the momentum correlations $(\tilde{g}^{(p)}(k_1 = k'_1, k_2 = k'_2, \dots))$ depend on the off-diagonals of the respective coordinate-space correlations $(g^{(p)}(x_1 \neq x'_1, x_2 \neq x'_2, \dots))$ and therefore yield additional information about the structure of the state that is not seen in the diagonals $g^{(p)}(x_1 = x'_1, x_2 = x'_2, \dots)$ [69]. In the following, we will drop the tilde for simplicity.

#### 4.3.1. 1st order momentum correlations

We normalize the 1st order RDM in momentum space with the local momentum densities and obtain the 1st order Glauber momentum correlation function

$$g^{(1)}(k'_1, k_1) = \frac{\rho^{(1)}(k_1, k'_1)}{\sqrt{\rho(k_1) \rho(k'_1)}} \quad (7)$$

Figure 3(a)–(d) shows $g^{(1)}(k'_1, k_1)$ for the various emerging phases. At $g_d = 0$, the system is fully coherent and $g^{(1)}$ has a constant value throughout k-space. This corroborates our observations of the coordinate-space correlations in figures 2(a), (e). As the interaction increases, $g^{(1)}(k'_1, k_1)$ shows an intricate diagonal pattern of alternating weak ($|g^{(1)}| \approx 1$) and strong ($|g^{(1)}| \approx 0$) correlations (figure 3(b)). With stronger interaction ($g_d = 1$), larger non-coherent regions form black stripes in figure 3(c). As the crystal phase is approached at
$g_d = 15$, the diagonal non-coherent regions broaden, forming an alternating high-low pattern (figure 3(d)) but with a larger spacing than figure 3(b). The stripe pattern formation is the result of the interaction between the BEC fragments. The spacing between the stripes relates inversely to the spatial distance between the individual fragments. In the MI state, the fragments are localized at each lattice site and their spatial distance matches the lattice constant. As the interaction increases, the bosons fragment further and they are no longer necessarily aligned to the lattice sites. In the CS, the effective lattice constant with respect to the fragments/bosons is roughly halved resulting in roughly doubled spacing between the stripes.

### 4.3.2. 2nd order momentum correlations

In figure 3(e)–(h) we show the 2nd order momentum correlation that offers insight into the 2nd order coherence of the state at distinct phases. The SF regime $g_d = 0$ (figure 3(e)) shows 2nd order coherence everywhere in $k$-space. At interaction strength $g_d = 0.1$ in the MI regime, a periodic bunching-anti-bunching pattern emerges (figure 3(f)). For stronger interactions the stripes in $g^{(2)}(k_1, k_2)$ become wider and lose contrast (figures 3(g)–(h)).

The form of the $k$-space correlations admits a simple interpretation. For weak interactions the gas is SF and this necessitates delocalization of the particles with complete $N$th-order coherence. This is reflected in the constant value of $|g^{(2)}| = 1$ that does not correlate or anti-correlate any pair-measurement of momenta $k_1, k_2$. At the MI state, finding one particle with momentum $k_1$ indicates the (almost) certainty to find particle 2 with $k_2 = k_1 + n\delta_k^{(\text{MI})}$, *n* integer. Well after the emergence of the TG-gas-alike state at $g_d \approx 1$, the stripe-like pattern re-emerges (see figure 3(h)), which can be interpreted as a fermionic diffraction pattern. Indeed, the existence of periodic patterns in $g^{(2)}(k_1, k_2)$ has been attributed to a phenomenon similar to Friedel oscillations in metals [27] where the (effective) fermion pair wavefunctions interfere. An examination of the structure of $g^{(2)}(k_1, k_2)$ reveals that the spacing of the stripes is $\delta_k^{(\text{MI})} \approx 2.3$ corresponding to a distance of $\delta_x^{(\text{MI})} \approx 1.36$ in real space, which is similar to the extension of each fermionized wavefunction, i.e. half a lattice-site width. This periodic pattern is more clearly reflected in the correlations at very large $g_d$ (figure 3(h)) in the crystal phase: there, the momentum distributions are expected—at least for the case of a harmonic confinement—to match the ones of interacting fermions [31]. When one particle is found at momentum $k_1$, the second one has a momentum of $k_2 = k_1 + m\delta_k^{(\text{CS})}$, with *m* integer. For our setup we find $\delta_k^{(\text{CS})} \approx 2\delta_k^{(\text{MI})}$. We stress that $\delta_k^{(\text{CS})}$ is a function of the interparticle interaction strength $g_d$ and does—unlike $\delta_k^{(\text{MI})}$—not relate to the underlying lattice spacing. In our setup, we find $\delta_k^{(\text{CS})} \approx 2\delta_k^{(\text{MI})}$, because of the hard-wall boundary conditions that we consider. Interestingly, the stripe-like periodic structure of $g^{(1)}$ and $g^{(2)}$ seen in the strongly interacting—fermionized—limit has been found for interacting fermions as well [94].

## 5. Effect of the lattice's depth

The two distinct phase transitions studied above show a very different dependence on the lattice depth $V$. The SF $\rightarrow$ MI transition shows a very strong dependence in $V$ while the MI $\rightarrow$ CS one is practically independent of $V$. In this section, we explore the dependence of the many-body correlations as well as the natural occupations on the lattice depth $V$ and assess how it relates fundamentally to the control of correlations in the system.

Figure 4(a) shows the first-order Glauber correlation function $|g^{(1)}(x, x')|$ for different lattice depths $V$. At $g_d = 0.02$ which is at the SF phase near the SF $\rightarrow$ MI transition, we see a stark dependence on $V$. For a small depth $V = 2$, the system is plainly SF and coherent. As the barrier is raised the ODLRO strongly reduces as can be seen from the decrease of the off-diagonal values of $|g^{(1)}|$ at $V = 6$. At $V = 10$ off-diagonal coherence is completely absent and a purely diagonal block structure emerges, indicating an MI phase. The MI phase persists with increasing $V$; a transition from MI to the CS is not seen (cf. plot for $V = 30$). Thus, keeping $g_d$ fixed to a small value at the SF phase, one can induce the transition from the coherent SF phase to a non-coherent MI phase solely by changing the lattice depth $V$.

In contrast, a transition from the MI to the CS or vice-versa cannot be induced by changing $V$; the MI $\rightarrow$ CS transition shows only a minimal $V$-dependence that can be seen in the lower panel of figure 4(a) that shows $|g^{(1)}|$ as a function of $V$ for $g_d = 15.0$. Here the system is in the crystal state and $|g^{(1)}|$ shows negligible change as $V$ is changed. The inter-site and intra-site correlations are thus minimally affected by the lattice depth as the system reaches the crystal phase.

This conclusion is supported by analyzing the effect of the lattice depth on the natural occupations. Figure 4(b) shows the natural occupation $\lambda_i$ as a function of interaction strength $g_d$ for different lattice depths $V$. For small interactions, there is a strong dependence on $V$. The larger the $V$, the smaller the critical $g_d$ at which the SF $\rightarrow$ MI transition occurs, which is revealed by the equal population of three natural orbitals ($\lambda_{1,2,3} \approx \frac{1}{3}$). The CS formation, indicated by an equal population of six orbitals, shows a minimal dependence on $V$.

These observations demonstrate the fundamental difference between the two distinct phase transitions as previously discussed. The SF $\rightarrow$ MI is an 'extrinsic' transition since it can be manipulated using the external
<figure><img src="image_7.png" /></figure>

(one-body) potential. The MI → CS one is, on the other hand, a purely ‘intrinsic’ transition, since it can be
induced solely by the two-body dipolar interactions and cannot be triggered via the one-body potential. As a
consequence, control of intersite correlations can be achieved through the manipulation of the external
potential while the intrasite correlations can be controlled primarily through the dipole–dipole interaction strength.

## 6. Conclusions

In this work, we have explored the many-body correlations of interacting ultracold dipolar bosons in optical lattices. As the dipolar interaction strength increases, the bosons transition from a superfluid to a Mott-insulating and, eventually, for stronger dipolar interactions, to a crystal phase. The phases are characterized by the natural orbitals and their populations. In the superfluid state, only one orbital is populated. In the Mott-insulating state, the orbitals that are equally populated are as many as the sites in the lattice. Finally, for even stronger dipolar interactions, the crystal phase emerges and the equally populated orbitals are as many as there are particles in the system. We thus verify a fundamental connection between the natural occupations and the strength of the dipole–dipole interaction [33]. Moreover, we examined the dependence of the natural orbitals on the interaction strength and showed how the densities in distinct quantum phases are built from them. Note that up to the value $g_d \approx 0.1$, all occupied natural orbitals topologically resemble the ones of the non-interacting system. Past $g_d \approx 0.1$, the natural orbitals become topologically distinct: they are torn or glued and not continuously deformed, i.e., the number of nodes in them changes as the interaction strength is increased.

As the interactions increase, dipolar bosons in optical lattices cover the full range of correlation properties: while the weakly interacting SF is entirely uncorrelated and the MI correlations are structured via the external lattice potential, the long-range interaction-dominated crystal phase represents a strongly correlated state of a many-body system structured by the interparticle interactions and not the lattice potential. We show how the normalized Glauber correlation functions of the ground state undergo characteristic changes as the interaction strength grows and the system transitions across different phases.

We comment on a fundamental difference between the SF → MI and the MI → CS transitions, which can be exploited for a control and management of the system’s correlations. The former transition is extrinsic, implying that it can be induced via the one-body potential. The latter is intrinsic, i.e. it is only accessed via the two-body interaction potential. Hence, precise control of the intersite (i.e. across different sites) correlations is possible through the manipulation of the depth of the lattice potential while the intrasite correlations can be controlled primarily by the dipolar interaction strength.

Recent experimental progress has propelled the development of a variety of methods for the measurement of spatial and momentum correlation functions thereby providing a detailed probe into the many-body properties of cold atom systems. In particular, time-of-flight measurements yield substantial information on correlation functions up to various orders [95–100] complemented by additional techniques such as noise pattern analysis [99], momentum microscopy [100], fluorescence imaging [96] etc.

Our current investigation of systems with dipole–dipole interactions provides a launchpad for further studies into the fundamental aspects of many-body correlations. A straightforward extension would be the calculation of higher-order correlation functions and also the investigation of lattices with incommensurate filling. Incommensurate systems are indeed drastically different from their commensurate counterparts [101, 102] and possess a stronger sensitivity to the exact particle numbers. Thus significantly different ground state properties are expected to be seen. A more general study that includes contact interactions along with dipolar ones could reveal the interplay between short- and long-range interactions with adjustable strengths and promises to show new multiscale phases [103].

## Acknowledgments

BC gratefully acknowledges financial support from the Department of Science and Technology, Government of India under the DST Inspire Faculty fellowship. MCT acknowledges financial support by the São Paulo state research foundation (FAPESP) and also CePOF/USP. AUJL acknowledges financial support by the Austrian Science Foundation (FWF) under grant No. P 32033, and the Wiener Wissenschafts- und TechnologieFonds (WWTF) project No MA16-066 ('SEQUEX'). Computation time on the HPC2013 cluster of the IIT Kanpur and the Hazel Hen cluster of the HLRS in Stuttgart are gratefully acknowledged.

## Appendix A. Results for larger systems

The main text showed the results for a triple-well system $S = 3$ with six bosons $N = 6$. In this section, we show similar spatial and momentum correlation characteristics for both larger lattice size $S = 5, N = 10$ and larger filling factor $S = 3, N = 9$ demonstrating that the analysis of the main text can be generalized for larger systems.
<figure><img src="image_8.png" /></figure>

**Figure A1.** (a)–(d): The 1st order spatial correlation function $|g^{(1)}(x_1, x_1')|$ and (e)–(h): the 2nd order correlation function $|g^{(2)}(x_1, x_2)|$, for the distinct emergent phases for $S = 3, N = 9$. At $g_d = 0$, the bosons are coherent, reflected in the uniform distribution of $|g^{(1)}|$ and $|g^{(2)}|$. For the Mott insulator ($g_d = 0.1$), the loss of intersite coherence but presence of intrasite coherence is reflected by a strong off-diagonal reduction and coherence triple-block-diagonal structure corresponding to $S = 3$ sites. Larger values of $g_d$ lead to the reduction of intra-well coherence and at the crystal phase, $g_d = 15$, coherence blocks centers at the location of each boson forming a nine-fold block structure corresponding to $N = 9$ bosons.

<figure><img src="image_9.png" /></figure>

**Figure A2.** (a)–(d): 1st order momentum correlations $g^{(1)}(k'_1, k_1)$ and (e)–(h): 2nd order momentum correlations $g^{(2)}(k_1, k_2)$ for the distinct emergent phases for $S = 3, N = 9$. Complete momentum coherence is seen in superfluid state for $g_d = 0$. At the Mott-insulator state($g_d = 0.1$), clean diagonal pattern of stripes with alternating weak and strong correlation is seen. At the crystal state, $g_d = 15$, a similar stripe pattern but with a different spacing is seen.

## A.1. Larger filling factors

To show the generality of our results for higher filling factors, we show the 1st and 2nd order spatial and momentum correlation functions for different emergent phases for $S = 3, N = 9$.

Figure A1 displays the 1st order spatial correlation function $|g^{(1)}(x_1, x_1')|$ and the 2nd order correlation function $|g^{(2)}(x_1, x_2)|$, for the distinct emergent phases for $S = 3, N = 9$.
<figure><img src="image_10.png" /></figure>

**Figure A3.** (a)–(d): The 1st order spatial correlation function $|g^{(1)}(x_1, x'_1)|$ and (e)–(h): the 2nd order correlation function $|g^{(2)}(x_1, x'_2)|$, for the distinct emergent phases for $S = 5, N = 10$. $g_d = 0$ shows complete coherence with uniform $|g^{(1)}|$ and $|g^{(2)}|$. At $g_d = 0.1$ (MI), the presence of intrawell coherence and absence of interwell coherence results in a five-fold block coherence structure corresponding to $S = 5$ sites. At the crystal state, $g_d = 15$, the complete localization of bosons is revealed in a ten-fold block structure as the coherence blocks centers on $N = 10$ bosons.

<figure><img src="image_11.png" /></figure>

**Figure A4.** (a)–(d): 1st order momentum correlations and (e)–(h): 2nd order momentum correlations for the distinct emergent phases for $S = 5, N = 10$. The momentum correlation shows a transition of complete coherence at $g_d = 0$ to a diagonal stripe pattern at the MI phase $g_d = 0.1$. A similar stripe pattern with different spacing is seen for crystal state at $g_d = 15$.

Figure A2 displays the 1st order momentum correlations $g^{(1)}(k'_1, k_1)$ and the 2nd order momentum correlations $g^{(2)}(k'_1, k_2)$ for the distinct emergent phases for $S = 3, N = 9$.

Both the 1st and 2nd order spatial and momentum correlation function for $S = 3, N = 9$ show similar characteristics as discussed in the main text for $S = 3, N = 6$. It can thus be concluded that the analysis of the spatial and momentum correlation functions of the main text is also valid for systems with larger filling.
## A.2. Larger lattice size

To demonstrate the generality of our analysis in a triple well for larger lattice size, we show the 1st and 2nd order spatial and momentum correlation functions for $S = 5, N = 10$.

Figure A3 displays the 1st order spatial correlation function $|g^{(1)}(x_1, x'_1)|$ and the 2nd order correlation function $|g^{(2)}(x_1, x_2)|$, for the distinct emergent phases for $S = 5, N = 10$.

Figure A4 displays the 1st order momentum correlations $g^{(1)}(k'_1, k_1)$ and (e)-(h): 2nd order momentum correlations $g^{(2)}(k_1, k_2)$ for the distinct emergent phases for $S = 5, N = 10$.

Both the 1st and 2nd order spatial and momentum correlation function for $S = 5, N = 10$ show similar characteristics as discussed in the main text for $S = 3, N = 6$. It can thus be concluded that the analysis of the spatial and momentum correlation functions of the main text is also valid for larger lattice size.

# Appendix B. One-body RDM

The reduced density matrix is important not only for its spectral decomposition but also from the perspective of its integral kernel $\rho^{(p)}(x_1, .. x_p, x'_1, .. x'_p) \equiv \langle x_1 .. x_p | \hat{\rho}^{(p)} | x'_1 .. x'_p \rangle$. The latter relates to the correlation and coherence properties of the system and is used to construct the generalized $p$th order correlation functions [69]

$$g^{(p)}(x_1, .. x_p, x'_1, .. x'_p) = \frac{\rho^{(p)}(x_1, .. x_p | x'_1, .. x'_p)}{\sqrt{\prod_{i=1}^{p} \rho(x_i) \rho(x'_i)}}. \quad (B1)$$

The 1st order correlations ($p = 1$) and hence the one-body aspects of the system are contained in the one-body RDM $\rho^{(1)}(x, x')$. Its diagonal kernel $\rho(x) \equiv \rho^{(1)}(x, x)$ is real-valued and gives the one-particle density or equivalently the probability of finding a particle in position $x$ irrespective of all other $N - 1$ positions. In contrast, the off-diagonal part, $\rho^{(1)}(x, x' \neq x)$, is in general complex and thus not directly experimentally observable. Physically, it represents the overlap of a particle state at position $x$ and $x'$.

<figure><img src="image_12.png" /><figcaption>**Figure B1.** One-body reduced density matrix $\rho^{(1)}(x, x')$ for various interaction strengths. For $g_d = 0$, the system is in a coherent superfluid state showing a checkerboard pattern with an intensity maximum at the origin. At $g_d = 0.1$, the system is in the MI state. The intersite coherence is reduced and reflected in the lower density values at the off-diagonals. At $g_d = 1$, we are in the Tonks-Girardeau regime and the density maxima in each well begins to split. At $g_d = 15$, we are deep in crystal phase and the complete splitting of the density maxima shows intrasite decoherence.</figcaption></figure>
The off-diagonal kernel of $\rho^{(1)}$ relates to the coherence properties of the system. For an infinite homogeneous system, non-vanishing off-diagonals $\rho^{(1)}(x, x')$ as $|x - x'| \to \infty$ for all pairs $(x, x')$ imply ODLRO and hence indicate coherence [104]. For a finite spatially bounded system, there is no strict ODLRO. Instead, for finite systems coherence is established when the off-diagonal of $\rho^{(1)}(x, x')$ is simply non-vanishing (without further requirements).

Figure B1 shows the one-body RDM $\rho^{(1)}(x, x')$ for various interaction strengths. In the limit of vanishing interaction, $g_d \to 0$ the bosons are condensed, and the state is SF. The comparatively high kinetic and potential energies lead to a greater population in the middle well. $\rho^{(1)}(x, x')$ thus shows a checkerboard pattern with a high value corresponding to the middle well and smaller values in the outer ones. As interaction is introduced, the bosons undergo localization in each well, whose degree increases with $g_d$ and is maximum at the MI phase for $g_d = 0.1$. This localization reflects in $\rho^{(1)}(x, x')$ which shows three diagonal maxima corresponding to each well, meaning that the probability of finding a particle in $x$ and $x'$ is significant only within each of the lattice sites. As the repulsive interaction increases, the particles inside each site begin to separate and thus localize individually. At the crystal phase $g_d = 15$, the bosons show complete intrasite localization with $\rho^{(1)}(x, x')$ showing $N = 6$ separate diagonal maxima corresponding to each particle.

# Appendix C. Momentum space RDM

## C.1. One-particle distribution

As in the coordinate-space density, the single-particle momentum density matrix relates to the momentum coherence of the system. Physically the first order momentum density matrix represents the probability of finding a particle with momenta $k$ and $k'$.

In the uncorrelated limit, $g_d = 0$, the momentum distribution is concentrated around $k = 0$ showing a checkerboard pattern with secondary momentum coherence at the reciprocal lattice sites (figure C1). When the interaction is introduced the coherence centers stretch diagonally forming a striped diagonal pattern with the principal diagonal having a maximum value in the MI phase. This reflects the localization in real space and, hence, delocalization in momentum space.

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /><figcaption>**Figure C2.** The two-particle momentum distribution $\rho^{(2)}(k_1, k_2)$ for various interaction strengths. At $g_d = 0$ we obtain a checkerboard pattern with peaked origin. In the MI state ($g_d = 0.1$) there appear three peaks over a background cloud, which are connected to the threefold fragmented state. At the crystal state the momentum distributions of strongly interacting bosons are expected to become equal to those of interacting fermions [31].</figcaption></figure>

At the MI phase $g_d = 0.1$, $\rho^{(1)}(k, k')$ shows three primary diagonal stripes pertaining to the three lattice sites. The side-bands form from stretching and joining the smaller coherence center at the reciprocal sites. There are also fainter bands at the edges formed from the stretching of the off-diagonal corners. Increasing interaction results initially, in a reduction of the distant stripes as can be seen at $g_d = 1.0$. However, as the interaction is further increased, the distant bands reappear while the nearer ones diminish. One can understand this behavior from the fact that there emerge two competing tendencies. The stronger repulsive interaction favors the enhancement of momentum correlations in the vicinity of $k = 0$. However, as the interaction increases, the condensate localizes in real space, resulting in the spreading of the momentum correlations.

At the crystal phase, $g_d = 15$, we get three well-separated stripes, signifying long-range momentum correlations.

### C.2. Two-particle distribution

The diagonal of the two-body momentum density (or, two-particle momentum distributions) gives the probability of finding one particle with momentum $k_1$ and another with momentum $k_2$. For $g_d \approx 0$ the $\rho(k_1, k_2)$ is concentrated near $k_1 = k_2 = 0$ with small contributions at $(0, \pm\pi/2)$ (figure C2). In this regime, as discussed, only one orbital contributes and the momentum distribution essentially reflects the diffraction of this orbital on the lattice.

For increasing interaction ($g_d = 0.1$) the discrete pattern now changes to a smeared distribution in the MI phase (figure C2). The momentum distributions concentrate on the diagonal forming three lobes over the continuous diagonal. This is contrasted to the real-space two-body density that has a correlation hole at $x_1 = x_2$ (see figure 2 of [33]). At $g_d = 15$, the system is at the crystal state and the six-fold fragmentation leads more structure at the diagonal, with further expansion and smearing of the off-diagonal contribution.

## References

[1] Baranov M A 2008 *Phys. Rep.* **464** 71

[2] Lahaye T, Menotti C, Santos L, Lewenstein M and Pfau T 2009 *Rep. Prog. Phys.* **72** 126401
[3] Griesmaier A, Werner J, Hensler S, Stuhler J and Pfau T 2005 *Phys. Rev. Lett.* **94** 160401

[4] Beaufils Q, Chicireanu R, Zanon T, Laburthe-Tolra B, Maréchal E, Vernac L, Keller J-C and Gorceix O 2008 *Phys. Rev.* A **77** 061601

[5] Lu M, Burdick N Q, Youn S H and Lev B L 2011 *Phys. Rev. Lett.* **107** 190401

[6] Aikawa K, Frisch A, Mark M, Baier S, Rietzler A, Grimm R and Ferlaino F 2012 *Phys. Rev. Lett.* **108** 210401

[7] Ni K-K, Ospelkaus S, de Miranda M H G, Pe'er A, Neyenhuis B, Zirbel J J, Kotochigova S, Julienne P S, Jin D S and Ye J 2008 *Science* **322** 231

[8] Takekoshi T, Reichsöllner L, Schindewolf A, Hutson J M, Le Sueur C R, Dulieu O, Ferlaino F, Grimm R and Nägerl H C 2014 *Phys. Rev. Lett.* **113** 205301

[9] Yi S and You L 2001 *Phys. Rev.* A **63** 053607

[10] Santos L, Shlyapnikov G V, Zoller P and Lewenstein M 2000 *Phys. Rev. Lett.* **85** 1791

[11] Góral K, Rzążewski K and Pfau T 2000 *Phys. Rev.* A **61** 051601

[12] Eberlein C, Giovanazzi S and O' Dell D H J 2005 *Phys. Rev.* A **71** 033618

[13] Góral K, Gajda M and Rzążewski K 2002 *Phys. Rev.* A **66** 023613

[14] Koch T, Lahaye T, Metz J, Fröhlich B, Griesmaier A and Pfau T 2008 *Nat. Phys.* **4** 218

[15] Ferrier-Barbut I, Kadau H, Schmitt M, Wenzel M and Pfau T 2016 *Phys. Rev. Lett.* **116** 215301

[16] Chomaz L, Baier S, Petter D, Mark M J, Wächtler F, Santos L and Ferlaino F 2016 *Phys. Rev.* X **6** 041039

[17] Wächtler F and Santos L 2016 *Phys. Rev.* A **94** 043618

[18] Macia A, Sánchez-Baena J, Boronat J and Mazzanti F 2016 *Phys. Rev. Lett.* **117** 205301

[19] Bruun G and Taylor E 2008 *Phys. Rev. Lett.* **101** 245301

[20] Cooper N and Shlyapnikov G 2009 *Phys. Rev. Lett.* **103** 155302

[21] Ni K, Ospelkaus S, Wang D, Quéméner G, Neyenhuis B, de Miranda M, Bohn J, Ye J and Jin D 2010 *Nature* **464** 1324

[22] Arkhipov A S, Astrakharchik G E, Belikov A V and Lozovik Y E 2005 *JETP Lett.* **82** 39

[23] Citro R, Orignac E, De Palo S and Chiofalo M 2007 *Phys. Rev.* A **75** 51602(R)

[24] De Palo S, Orignac E, Citro R and Chiofalo M 2008 *Phys. Rev.* B **77** 212101

[25] Pedri P, De Palo S, Orignac E, Citro R and Chiofalo M 2008 *Phys. Rev.* A **77** 15601

[26] Zöllner S, Brunn G M, Pethick C J and Reimann S M 2011 *Phys. Rev. Lett.* **107** 035301

[27] Zöllner S 2011 *Phys. Rev.* A **84** 063619

[28] Maik M, Buonsante P, Vezzani A and Zakrzewski J 2011 *Phys. Rev.* A **84** 053615

[29] Astrakharchik G E, Morigi G E, De Chiara G and Boronat J 2008 *Phys. Rev.* A **78** 063622

[30] Astrakharchik G E and Lozovik Y E 2008 *Phys. Rev.* A **77** 013404

[31] Deuretzbacher F, Cremon J C and Reimann S M 2010 *Phys. Rev.* A **81** 063616

[32] Chatterjee B, Brouzos I, Cao L and Schmelcher P 2013 *J. Phys. B: At. Mol. Opt. Phys.* **46** 085304

[33] Chatterjee B and Lode A U J 2018 *Phys. Rev.* A **98** 053624

[34] Góral K, Santos L and Lewenstein M 2002 *Phys. Rev. Lett.* **88** 170406

[35] Dalla Torre E G, Berg E and Altman E 2006 *Phys. Rev. Lett.* **97** 260401

[36] Deng X and Santos L 2011 *Phys. Rev.* B **84** 085138

[37] Menotti C, Trefzger C and Lewenstein M 2007 *Phys. Rev. Lett.* **98** 235301

[38] Capogrosso-Sansone B, Trefzger C, Lewenstein M, Zoller P and Pupillo G 2010 *Phys. Rev. Lett.* **104** 125301

[39] Capello M, Becca F, Fabrizio M and Sorella S 2007 *Phys. Rev. Lett.* **99** 056402

[40] Petrov D S, Shlyapnikov G V and Walraven J T M 2000 *Phys. Rev. Lett.* **85** 3745

[41] Dunjko V, Lorent V and Olshanii M 2001 *Phys. Rev. Lett.* **86** 5413

[42] Paredes B, Widera A, Murg V, Mandel O, Fölling S, Cirac I, Shlyapnikov G V, Hänsch T W and Bloch I 2004 *Nature* **429** 277

[43] Deuretzbacher F, Bongs K, Sengstock K and Pfannkuche D 2007 *Phys. Rev.* A **75** 013614

[44] Jaksch D, Bruder C, Cirac J I, Gardiner C W and Zoller P 1998 *Phys. Rev. Lett.* **81** 3108

[45] Greiner M, Mandel O, Esslinger T, Hänsch T W and Bloch I 2002 *Nature* **415** 39

[46] Girardeau M 1960 *J. Math. Phys.* **1** 516

[47] Serwane F, Zürn G, Lompe T, Ottenstein T B, Wenz A N and Jochim S 2011 *Science* **332** 336

[48] Murmann S, Bergschneider A, Klinkhamer V M, Zürn G, Lompe T and Jochim S 2015 *Phys. Rev. Lett.* **114** 080402

[49] Murmann S, Deuretzbacher F, Zürn G, Bjerlin J, Reimann S M, Santos L, Lompe T and Jochim S 2015 *Phys. Rev. Lett.* **115** 215301

[50] Peter D, Pawlowski K, Pfau T and Rzążewski K 2012 *J. Phys.* B **45** 225302

[51] Lahaye T, Pfau T and Santos L 2010 *Phys. Rev. Lett.* **104** 170404

[52] Dell'Anna L, Mazzarella G, Penna V and Salasnich L 2013 *Phys. Rev.* A **87** 053620

[53] Xiong B and Fischer U R 2013 *Phys. Rev.* A **88** 063608

[54] Gallemi A, Guilleumas M, Mayol R and Sanpera A 2013 *Phys. Rev.* A **88** 063645

[55] Gallemi A, Queralto G, Guilleumas M, Mayol R and Sanpera A 2016 *Phys. Rev.* A **94** 063626

[56] Biedroń K, Łacki M and Zakrzewski J 2018 *Phys. Rev.* B **97** 245102

[57] Alon O E, Streltsov A I and Cederbaum L S 2008 *Phys. Rev.* A **77** 033613

[58] Cao L, Mistakidis S I, Deng X and Schmelcher P 2017 *Chem. Phys.* **482** 303

[59] Mistakidis S I, Cao L and Schmelcher P 2015 *Phys. Rev.* A **91** 033611

[60] Koutentakis G M, Mistakidis S I and Schmelcher P 2017 *Phys. Rev.* A **95** 013617

[61] Roy R, Gammal A, Tsatsos M C, Chatterjee B, Chakrabarti B and Lode A U J 2018 *Phys. Rev.* A **97** 043625

[62] Bera S, Chakrabarti B, Gammal A, Tsatsos M C, Lekala M L, Chatterjee B, Lévèque C and Lode A U J 2018 arXiv:1806.02539

[63] Lode A U J, Tsatsos M C, Fasshauer E, Weiner S E, Lin R, Papariello L, Molignini P and Lévèque C 2018 MCTDH-X: The time-dependent multiconfigurational Hartree for indistinguishable particles software (http://ultracold.org)

[64] Fasshauer E and Lode A U J 2016 *Phys. Rev.* A **93** 033635

[65] Lode A U J 2016 *Phys. Rev.* A **93** 063601

[66] Lühmann D S, Bongs K, Sengstock K and Pfannkuche D 2008 *Phys. Rev.* A **77** 023620

[67] Sinha S and Santos L 2007 *Phys. Rev. Lett.* **99** 140406

[68] Cai Y, Rosenkranz M, Lei Z and Bao W 2010 *Phys. Rev.* A **82** 043623

[69] Sakmann K, Streltsov A I, Alon O E and Cederbaum L S 2008 *Phys. Rev.* A **78** 023615

[70] Glauber R J 1963 *Phys. Rev.* **130** 2529

[71] Penrose O and Onsager L 1956 *Phys. Rev.* **104** 576

[72] Nozières P and Saint James D 1982 *J. Phys.* **43** 1133
[73] Spekkens R W and Sipe J E 1999 Phys. Rev. A 59 3868

[74] Mueller E J, Ho T L, Ueda M and Baym G 2006 Phys. Rev. A 74 033612

[75] Alon O E, Streltsov A I and Cederbaum L S 2005 Phys. Rev. Lett. 95 030405

[76] Fischer U R, Lode A U J and Chatterjee B 2015 Phys. Rev. A 91 063621

[77] Bloch I, Dalibard J and Zwerger W 2008 Rev. Mod. Phys. 80 885

[78] Boéris G et al 2016 Phys. Rev. A 93 011601(R)

[79] Astrakharchik G E, Krutitsky K V, Lewenstein M and Mazzanti F 2016 Phys. Rev. A 93 021605(R)

[80] Meinert F, Knap M, Kirilov E, Jag-Lauber K, Zvonarev M B, Demler E and Nägerl H-C 2017 Science 356 945

[81] Rigol M and Muramatsu A 2004 Phys. Rev. A 70 031603(R)

[82] Girardeau M D and Astrakharchik G E 2012 Phys. Rev. Lett. 109 235305

[83] Sakmann K, Streltsov A I, Alon O E and Cederbaum L S 2011 New J. Phys. 13 043003

[84] Streltsov A I, Cederbaum L S and Moiseyev N 2004 Phys. Rev. A 70 053607

[85] Lode A U J and Bruder C 2017 Phys. Rev. Lett. 118 013603

[86] Krönke S and Schmelcher P 2015 Phys. Rev. A 92 023631

[87] Mistakidis S I, Katsimiga G C, Kevrekidis P G and Schmelcher P 2018 New J. Phys. 20 043052

[88] Cirac J I, Gardiner C W, Naraschewski M and Zoller P 1996 Phys. Rev. A 54 R3714

[89] Naraschewski M, Wallis H, Schenzle A, Cirac J I and Zoller P 1996 Phys. Rev. A 54 2185

[90] Javanainen J and Wilkens M 1997 Phys. Rev. Lett. 78 4675

[91] Naraschewski M and Glauber R J 1999 Phys. Rev. A 59 4595

[92] Yang C N 1962 Rev. Mod. Phys. 34 694

[93] Emary C, Pölzl C, Carmele A, Kabuss J, Knorr A and Brandes T 2012 Phys. Rev. B 85 165417

[94] Brandt B B, Yannouleas C and Landman U 2017 Phys. Rev. A 96 053632

[95] Bouchoule I, Arzamasovs M, Kherunstsyans K V and Gangardt D M 2012 Phys. Rev. A 86 033626

[96] Perrin A, Bucker R, Manz S, Betz T, Koller C, Plisson T, Schumm T and Schmiedmayer J 2012 Nat. Phys. 8 195

[97] Dall R G, Manning A G, Hodgman S S, RuGway W, Kherunstsyans K V and Truscott A G 2013 Nat. Phys. 9 341

[98] Langen T, Erne S, Geiger R, Bauer B, Schweigler T, Kuhnert M, Rohinger W, Mazets I E, Gasenzer T and Schmiedmayer J 2015 Science 348 207

[99] Fang B, Johnson A, Roscilde T and Bouchoule I 2016 Phys. Rev. Lett. 116 050402

[100] Hodgman S S, Khakimov R I, Lewis-Swan R J, Truscott A G and Kherunstsyans K V 2017 Phys. Rev. Lett. 118 240402

[101] Brouzos I, Zöllner S and Schmelcher P 2010 Phys. Rev. A 81 053613

[102] Cinti F 2016 J. Low Temp. Phys. 182 513

[103] Pizzardo M, Mazzarella G and Salasnich L 2016 J. Low Temp. Phys. 185 59

[104] Yang C N 1962 Rev. Mod. Phys. 34 694

[105] Streltsov A I, Alon O E and Cederbaum L S 2007 Phys. Rev. Lett. 99 030402

[106] Lode A U J, Chakrabarti B and Kota V K B 2015 Phys. Rev. A 92 033622

[107] Weiner S E, Tsatsos M C, Cederbaum L S and Lode A U J 2017 Sci. Rep. 7 40122

[108] Tsatsos M C and Lode A U J 2015 J. Low Temp. Phys. 181 171

[109] Lode A U J, Klaiman S, Alon O E, Streltsov A I and Cederbaum L S 2014 Phys. Rev. A 89 053620

[110] Lode A U J and Bruder C 2016 Phys. Rev. A 94 013616

[111] Lode A U J et al 2018 New J. Phys. 20 055006