Strongly dipolar gases in a one-dimensional lattice: Bloch oscillations and matter-wave localization

Three-dimensional quantum gases of strongly dipolar atoms can undergo a crossover from a dilute gas to a dense macrodroplet, stabilized by quantum fluctuations. Adding a one-dimensional optical lattice creates a platform where quantum fluctuations are still unexplored, and a rich variety of new phases may be observable. We employ Bloch oscillations as an interferometric tool to assess the role quantum fluctuations play in an array of quasi-two-dimensional Bose-Einstein condensates. Long-lived oscillations are observed when the chemical potential is balanced between sites, in a region where a macrodroplet is extended over several lattice sites. Further, we observe a transition to a state that is localized to a single lattice plane$-$driven purely by interactions$-$marked by the disappearance of the interference pattern in the momentum distribution. To describe our observations, we develop a discrete one-dimensional extended Gross-Pitaevskii theory, including quantum fluctuations and a variational approach for the on-site wavefunction. This model is in quantitative agreement with the experiment, revealing the existence of single and multisite macrodroplets, and signatures of a two-dimensional bright soliton.

The dipole-dipole interaction (DDI) between magnetic atoms in an ultracold quantum gas has been key to the discovery of supersolids [1][2][3] and macrodroplets [4,5], new states of matter with extremely intriguing and counter-intuitive properties [6,7].Macrodroplets are macroscopic quantum states that behave in many ways like liquid droplets [4,5,8,9].They are at least an order of magnitude denser than normal Bose-Einstein condensates (BECs), and can be self-bound.They exist in a parameter regime in which mean-field theories predict the collapse of the entire system when the attractive dipolar interactions overcome the repulsive contact interactions.Instead, the system remains surprisingly stable thanks to the so-called quantum fluctuations, thus providing one of the rare examples where beyond-meanfield interactions substantially change the ground state of the system [10,11].Although the functional form of the beyond-mean-field term, otherwise known as the Lee-Huang-Yang (LHY) correction [12], is still subject to intense study and debate [13,14], its importance is now undoubted.Isolating beyond-mean-field effects may be crucial to settle disputes on its validity, particularly in dipole dominated systems; however, it is very difficult to have access to individual interaction contributions.Though, the differing atom number scaling between mean-field and LHY contributions provide a promising method to differentiate between them.
Optical lattices enable powerful interferometric approaches to, e.g., measure with high precision the zerocrossing of the scattering length or of the mean-field interaction with the so-called Bloch oscillation (BO) technique [15][16][17][18], and to achieve an accurate determination of the background scattering length via lattice spectroscopy in Hubbard models [4,19,20].Moreover, the presence of the lattice itself may change completely the phase diagram of the system, as shown in seminal experiments with contact interacting gases [21][22][23].Unique phenomena are predicted with the addition of long-range DDIs [24,25].Experiments with lattice-confined atomic dipolar gases have already shown important results, e.g., the realization of extended Bose-Hubbard models [26] and spin models [27][28][29][30] in three-dimensional (3D) lattices.In 2D lattices, forming quasi-1D tubes, suppression of dipolar relaxation [31] and the controlled breakdown of integrability [32] have been observed.Instead, up to now, 1D lattices, forming an array of quasi-2D layers, have been used with large wavelengths to load a single pancake trap [33], or multi-layer traps to study the role of DDI in the stability against collapse [34].Further, theoretical proposals have suggested that the DDI between layers not only can lead to modifications within each layer [35][36][37][38] but also to inter-layer bound states [39][40][41].Other works predict the existence of bright-soliton structures along the lattice [42] or anisotropic on-site solitons [43,44].However, those proposals lack the important stabilization mechanism given by the LHY term, which is known to provide many new phases in continuous systems (e.g.harmonically trapped), opening up many questions: What is the ground state of an attractive dipolar gas in a 1D lattice potential?Can droplets be delocalized over many lattice planes?Will solitonic solutions continue to exist?
In the present work, we study an erbium dipolar gas in a 1D optical lattice with dominantly attractive DDI.We employ BOs as an interferometric tool to probe the interaction contributions of the system, and to isolate the role of beyond-mean-field effects.We find long-lived oscillations, associated with a minimum in the dephasing rate, close to the cancellation point between meanfield and beyond-mean-field interactions, and at scatter- ing lengths significantly shifted from the expected meanfield result.We develop a discrete effective 1D extended Gross-Pitaevskii equation (eGPE) with variational transverse widths [45,46].We find that this minimum occurs when the chemical potentials on each site are equal, not the energies-as has been employed successfully in contact interaction dominated systems [17,18]-due to the difference in density scaling between the interactions.The close correspondence between theory and experiment shows the validity of the LHY prediction, even while highly inhomogeneous densities are expected to break the local density approximation [12].Moreover, we see that for low scattering lengths the system undergoes a structural transition to a single localized 2D plane, signifying an important new way to generate systems in reduced geometries through varying the interactions alone.Finally, using our theoretical model we produce a full phase diagram of the system, revealing the impact of the LHY contribution to the predicted 2D anisotropic soliton state [43], which is instead morphed into a droplet solution at high atom numbers.Though, promisingly, we still find soliton-like solutions exist.
In the experiment, we prepare a degenerate dipolar gas of erbium atoms in a one-dimensional optical lattice as follows.We start with a dipolar quantum gas of 5 × 10 4 spin-polarized 166 Er atoms confined in a cigar-shaped optical dipole trap [47] elongated along y.Typical BEC fractions range from 60% to 80%.The dipolar length for 166 Er is fixed at a dd = 66.5 a 0 , where a 0 is the Bohr radius.We tune the contact interaction between atoms and therefore the s-wave scattering length, a s , via Feshbach resonances [4,[48][49][50] by changing the absolute value of a bias magnetic field |B|.We fix the orientation of B to be along the weak axis (y) of the trap, making the DDI dominantly attractive [4,7].
Once the harmonically-trapped cloud is prepared at the desired a s , we switch on a 1D optical lattice, aligned along the gravity direction (z); see Fig. 1(a).The vertical lattice is created by retro-reflecting a λ = 1064 nm laser beam.We load the planes by exponentially increasing the lattice depth V 0 to 8 E rec in 20 ms, where E rec = 2 k 2 /2m = h × 10.5 kHz.Here, = h/2π is the reduced Planck's constant (h), m is the mass of 166 Er atoms and k = 2π/λ is the wave-vector of the lattice.The 1D lattice forms an array of tightly confined quasi-2D planes with a trap frequency along the tight direction ω z 2π × 6 kHz, corresponding to an harmonic oscillator length z ho = 100 nm.The tunnelling rate, J, between planes is about h × 33 Hz.For these 1D lattice parameters, ω z > k B T and the system is kinematically 2D [51].
We first aim at inducing Bloch oscillations to interferometrically assess the role of beyond-mean-field effects and test the validity of the LHY term.We thus suddenly switch off the dipole trap and let the system evolve in the combined lattice and gravitational potential for a variable hold time t h .Finally, using standard absorption imaging after 30 ms of time-of-flight (TOF), we record the evolution of the momentum distribution and extract the position of the main peak, q max , as a function of t h .Figure 1(b) shows an exemplary set of absorption images during a single Bloch period T BO .We observe the key paradigm of BOs, i.e. the linear increase of the mean momentum due to the acceleration and the Bragg reflection occurring at the border of the Brillouin zone [52], well described by fitting a sawtooth function to q max .
The high sensitivity of BOs to interactions [17,53] clearly appears by tracing the evolution for two different a s (see Fig. 1(c,d)), as the interaction dependence is encoded into the dephasing rate.For a contact-dominated gas (a dd < a s = 90 a 0 , Fig. 1(c)), we see that the BOs vanish within a few T BO .On the contrary, decreasing a s , and thereby going into the regime where contact interactions and DDI nearly compensate each other (a s = 60 a 0 , Fig. 1(d)), we observe persisting oscillations for more than 25 Bloch cycles, set by our limited observation time [50].To systematically study this effect, we repeat the BO measurements for different values of a s , and extract the corresponding dephasing rate γ [50].As shown in Fig. 2(a), we observe a resonant-type behavior with γ showing a pronounced dip with a minimum at a s = 61 a 0 .This minimum is clearly different to the point a s ≈ a dd , where the variance of the mean-field energies across different lattice sites cancel [18], which would be expected from previous observations [17,53].
To get further insight on the origin of the minimum, we develop a discrete effective 1D eGPE, inspired by the close correspondence between predictions from discrete models and experimental observations in non-dipolar [54,55] and weakly dipolar [18] BECs.We separate the 3D wavefunction into radial and axial contributions, allowing for a variational anisotropic radial width and thus maintaining the 3D character [45].Along the lattice direction (z), we further decompose the wavefunction, ψ(z, t), as a sum of Wannier functions w(z) of the lowest energy band over all lattice sites: ψ(z, t) = √ N j c j (t)w(z − z j ), where N is the atom number and c j (t) the complex wavefunction amplitude on lattice site j, leading to a set of discrete effective 1D eGPEs, each including mean-field and beyond mean-field interactions.For the beyond-mean-field interaction, the 3D form of the LHY still fully applies since the contact interaction energy exceed the confinement energy scale [50,56].However, our system may also open to further studies on the 2D to 3D crossover of the LHY.We solve these equations coupled to a minimization of the energy functional with respect to the variational parameters to determine the ground states, benchmarking them against the full 3D theory.We then perform dynamic simulations of the expected time evolution [50], giving an accurate dephasing rate (solid line) in Fig. 2(a) without free parameters.
In previous studies, the point of minimum dephasing was found to occur when the mean-field interaction energies vanish or cancel.We isolate the mean-field contribution by removing beyond-mean-field effects from our simulations (dashed line in Fig. 2(a)), predicting a minimum at a s ≈ a dd .However, this is in clear contradiction with our experimental observations by a shift of 6a 0 and a different overall shape due to the different scaling of the LHY term with the density.Without LHY, the cancellation of mean-field energies, E j MF , is equivalent to the cancellation of onsite chemical potentials, given by µ j = 2E j MF /|c j | 2 .Note, µ j dictates the wavefunction phase winding on each site through c j = |c j |e −iµj t/ .Reintroducing quantum fluctuations, we obtain , where the 5/2 appears due to the |c j | 5 density scaling in the beyondmean-field energy (E j BMF ). Figure 2(b) shows µ j from the ground state calculation for four scattering lengths, additionally indicating the contribution of the LHY correction.
We observe that the point of minimal dephasing in the experiment is close to the point where the variance of µ j is minimized [57].Indeed, within a semi-analytic approximation (see Ref. [50] for details), we find a direct relationship between γ and µ j , which reads γ ∝ |µ 1 − µ 0 | when 3 lattice sites (j = −1, 0, 1) are occupied.This model can be extended to 5 lattice sites, giving the dotdashed line (Fig. 2(a)) which reproduces very well the system behaviour [50].Interestingly, measuring the dephasing rate through the chemical potential is ubiquitous to systems with arbitrary interaction potentials.Surprisingly, by further decreasing the scattering length below 57 a 0 , no BOs nor interference peaks are visible anymore.We observe at the initial instant (t h = 0 ms) that the momentum distribution is already spread over the entire first Brillouin zone.To quantify this, we study the contrast, C, of the interference pattern of the initial momentum distribution as a function of a s , see Figure 3(a).We extract C, defined as the amplitude of the momentum peaks at ±2 k relative to the zero momentum peak, from the Fourier analysis of the TOF images [50].For large a s , we observe the typical matterwave interference pattern, as expected from a coherent state populating several lattice planes (see inset) [55].As we lower a s , C first remains fairly constant.For a s below a certain critical value a * s ≈ 57 a 0 , we observe a sudden loss of the interference pattern with a sharp decrease of C to almost zero.
Remarkably, we observe that this interaction-driven process is reversible.To test the restoring of the interference pattern, we employ the following protocol [50]: In brief, we first prepare the system in the lattice at constant and large a s (a s = 69(2) a 0 ).We then ramp down a s below a * s (a s = 56(2) a 0 ) in 20 ms and wait until C stabilizes to a small value; see Fig. 3(b).Note that the interference pattern disappears after about 10 ms, which is on the order of the tunneling time h/J between two neighboring lattice sites.At this point, we quench a s back to its initial value and probe the time evolution of the system towards its new equilibrium state.On a similar timescale, we observe the reappearance of the interference pattern with an increase of C, which then saturates to about 60% of its initial value [58].For comparison, we also show the data without inverting the field ramp.
The observed broad distribution in reciprocal space suggests that the system ground state has undergone a structural change, with the macroscopic wavefunction localized in one lattice plane.To verify this interpretation, we calculate the ground state of the system as a function of a s .When the repulsive contact interaction dominates (a s > a dd ), we find an array of BECs occupying approximately three to five lattice planes; see insets Fig. 3(a).In contrast, when the relative strength of the attractive dipolar interaction with respect to the other terms in the Hamiltonian is increased, the system reaches a critical point.Here, it undergoes a phase transition to a quasi-2D state, in which all atoms are localized into a single lattice plane to minimize their energy.This purely interactiondriven phase transition-somewhat reminiscent of a continuous version of a superfluid to Mott insulator transition [59]-is stabilized by quantum fluctuations (LHY), preventing the subsequent collapse of the system [42,60].The predicted critical point occurs exactly where we observe the disappearance of the interference pattern in the experiments.We find an overall excellent agreement between the measured and the calculated C from both the discrete 1D model and the 3D theory without any free fitting parameters, except for a rescaling factor to the contrast amplitude to account for the thermal atoms in the experiment.
The observation of this phase transition to a quasi-2D localized state driven by interactions points to the existence of a rich variety of phases.The importance of the LHY correction and its peculiar density scaling motivate us to investigate the properties of the ground state as a function of a s and atom number to identify distinct phases in this unique setting.For this, we employ our discrete model to derive a full phase diagram; see Fig. 4(a).To investigate the boundness of the states, we assess the impact of the radial harmonic trap on the minimum of the variational energy, which is a function of the radial widths l x and l y .At large scattering lengths, as expected, we find a stable delocalized BEC phase, where the total interaction energy (mean-field + LHY) is positive.The state is trap-bound, meaning that there is no energy minimum without the radial harmonic confinement; inset of Fig. 4(a).
Reducing a s , we find an energy minimum even without the radial harmonic trap (colored region in the phase diagram).These quasi-2D self-bound solutions (the lattice still provides axial confinement) are either extended over several sites (lighter color) or localized to a single plane (darker color).In the literature, there are two paradigmatic examples of self-bound objects with attractive mean-field energy: droplets and solitons.Droplets can exist in one, two or three dimensions and are stabilized through the LHY correction [7].Stable bright solitons only exist in quasi-1D systems with attractive contact interactions and are stabilized against collapse purely by kinetic energy.In the search for solitons in higher dimensions, theoretical studies have suggested that the DDI could stabilize such 2D solutions [43,44].To the best of our knowledge, there have been no studies on the effect the LHY correction has on this prediction, nor experimental observation.In the present case, where many interactions and kinetic energy compete, a classification of self-bound solutions is much less straightforward.As a crucial distinction between a soliton and a droplet, we use the scaling of the system width with atom number.The soliton width (along the collapse direction) scales inversely with increasing atom number [61], while in contrast, the droplet size increases in all directions with N [62], as predicted in a quasi-1D setting [63].We use this distinction to draw a boundary between the two phases, observing a phase transition at around 5000 atoms, for both single-site and multi-site solitons.The overlaying of our measurements (Fig. 3(a)) onto the phase diagram suggests that the experiments have already reached the interesting regimes of both 2D self-bound droplet and dipolar solitons.This opens the door to future experimental investigation on the self-bound nature and properties of these new 2D phases.
In conclusion, we theoretically and experimentally investigate the behavior of a strongly dipolar quantum gas in a 1D optical lattice.We employ BOs and characterize their dephasing rate as a function of a s .We observe a minimum in the dephasing shifted 6 a 0 away from the purely mean-field prediction, providing an interferometric measure of the beyond-mean-field contribution.For low enough a s , the system enters into a quasi-2D state which is localized onto a single lattice plane, providing a genuine interaction-driven path to reach reduced dimensions in dipolar gases.Using our developed discrete theory model, we derive a full phase diagram which confirms the observed localization transition.This also reveals signatures of quasi-2D self-bound dipolar droplet solutions, and the long sought-after 2D anisotropic dipolar soliton, first predicted in Ref. [43] (see also [44,64]).Our work paves the way for future studies of the soliton-to-droplet crossover in a dipolar gas, as observed in a Bose-Bose gas [65], and of the "solitonic" nature [66] of dipolar solitary waves [67][68][69][70][71].In this work, we use an extended Gross-Pitaevskii theory for direct comparison to our experimental results.We employ both the standard three-dimensional form of the extended Gross-Pitaevskii equation (eGPE) and derive a discrete effective one-dimensional eGPE.Starting with the three-dimensional case, our system can be described by the 3D eGPE of the form [4,[72][73][74] where the wavefunction Ψ is normalized to the total atom number N = d 3 x |Ψ| 2 .The atoms are confined in a harmonicpotential V harm = ξ=x,y,z 1 2 mω 2 ξ ξ 2 with single particle mass m and trap frequencies ω ξ , together with the lattice potential V latt = sE rec sin 2 (kz) where s is the tunable lattice depth in multiples of the recoil energy E rec and k = 2π/λ is the lattice spacing in reciprocal space.The mean-field interaction contributions are g = 4π 2 a s /m for the contact interaction, governed by the s-wave scattering length a s , and the long-ranged anisotropic dipolar interaction potential U dd ( x) = 3 2 a dd /m 1 − 3 cos 2 θ /| x| 3 , where a dd = µ 0 µ 2 m m/12π 2 with magnetic moment µ m and θ is the angle between the polarization axis (y-axis) and the vector between neighboring atoms.We also include beyondmean-field effects through the quantum fluctuations term [12], which depends on the relative strength between the dipolar and short-ranged interactions ε dd = a dd /a s .Finally, F ext = g grav m de-notes the external force exerted on the system by gravity.
In this work, we employ the imaginary time-evolution technique on Eq. (S1) in order to find stationary solutions for the wavefunction in the lattice, without gravity.For various atom numbers and scattering lengths, we use a numerical grid of lengths (L x , L y , L z ) = (6, 33.3, 6) µm, with corresponding grid points 128 × 256 × 128.The dipolar term is efficiently calculated in momentum space, and we use a cylindrical cut-off in order to negate the effects of aliasing from the Fourier transforms [75].
To derive the effective one-dimensional model, we follow Ref. [45] by assuming a wavefunction decomposition (S2) with variational parameters l and η representing the width of the radial wavefunction and the anisotropy of the state, respectively.Integrating out the transverse directions (x, y) in Eq. (S1) upon substitution of the ansatz above gives the continuous quasi-one-dimensional eGPE, which when combined with a variational minimization of the energy functional to find (l, η) gives close agreement to the full 3D eGPE [45].We further decompose the longitudinal wave function ψ(z, t) into a sum of Wannier functions w(z) of the lowest energy band over all lattice sites for complex amplitudes c j , and positions of lattice minima z j = (λ/2)j.For deep enough lattices, the Wannier functions are well approximated by Gaussians of the form After multiplying on the left by c * j and integrating over z, we obtain a set of discrete effective one-dimensional eG-PEs with the reduced effective one-dimensional parameters γ 1D QF = 2 3/2 /(5π 3/2 l 2 l latt ) 3/2 γ QF and g 1D = g/((2π) 3/2 l 2 l latt ).Here, J denotes the tunneling rate between two neighboring lattice sites.The dipolar interaction coefficients between lattice sites j and k depend both on the separation |j − k|, and non-trivially on the size l and anisotropy η of the transverse cloud.For the variational minimization, we generate an interpolating function for a sensible range of (l, η) and separations up to |j − k| = 6 via where Ψ 0 ( x , l, η) = Φ(x, y, l, η)w(z) [see Eqs.(S2) and ( S3)].This allows us to simply look up the values of U dd |j−k| without having to recalculate for every time step during the energy minimization.We note that the energy contribution rapidly declines for separations larger than 2 sites, and find that 6 is more than sufficient to quantitatively describe the physics.
To find the stationary state solution of Eq. (S4) (without gravity) we employ an imaginary time-evolution in combination with an optimization scheme, aiming to find the state which minimizes the total energy functional where c = (c 1 , c 2 , . . ., c n ) for n total lattice sites.Here, E ⊥ [l, η] gives the energy contribution from the transverse variational wave function, which reads The latter term of Eq. (S6) gives the discrete energy functional for the amplitudes c j , which includes the tunneling and all interaction terms Starting from an initial distribution of the amplitudes c j we first determine the variational parameters (l, η), which is done via an optimization scheme minimizing Eq. (S8).Subsequently, we evolve the amplitudes in imaginary time using Eq.(S4) and repeat this process until we find the minimum of the total energy function Eq. (S6).
In Fig. S1 we assess the different interaction energy contributions to Eq. (S6) for a range of scattering lengths.For a s > a dd the total interaction energy is positive, and it corresponds to a dilute BEC.Following a s to smaller values all interaction contributions are almost constant, until at around a s = 60 a 0 there is a phase transition from the BEC to droplet state, as identified in Fig. 4 of the main text.This sharp gradient ceases at around a s = 55 a 0 , where the atoms are localized to a single lattice plane.Note that although the DDI offsite energy is typically only 10% of the onsite counterpart, it constitutes a significant contribution to the total interaction energy in the system, shifting the BEC to droplet crossover and localization transitions by a few a 0 .
Once we have the ground state of the system, we employ the discrete effective one-dimensional eGPE in realtime to simulate the Bloch oscillations in the presence of gravity.

2D TO 3D CROSSOVER
The dimensionality of the system is known to highly influence the size and even the sign of the beyond-meanfield contribution, in both Bose-Bose [76][77][78] and dipolar [56,79,80] gases.Here, we assess the validity of employing the full 3D LHY correction to our system.Following Ref. [56], we define the dimensionless parameter ξ = gn/ 0 -dependent on the contact interactions g, peak 3D density n, and the confinement energy scale 0 = 2 π 2 /2mz 2 ho -that indicates which dimensionality regime our system is in.If ξ 1 we are safe to use the 3D LHY term, whereas if ξ 1 the 2D solution deviates from the 3D one.Deep in the localized droplet regime, where the peak density is on the order of 10 22 m −3 , we find ξ ≈ 2, and the 3D LHY as used throughout this work is valid.Even at large scattering lengths, where the peak density is closer to 5 × 10 20 m −3 , we find ξ ≈ 0.5, which introduces an error of less than 5% between the 2D and 3D LHY terms [56].In this limit, the 2D LHY term may be more appropriate, however in the dilute BEC phase the impact of the LHY is minimal.

ANALYTIC MODEL OF DEPHASING
Starting from the discrete 1D eGPE we decompose the coefficients c j into amplitude and phase as c j = |c j | exp(−iφ j ), and then integrate Eq. (S4) in time to give with onsite chemical potentials µ j , and where we have also assumed that F ext d J such that the amplitudes |c j | are frozen.
Following Ref. [81], we write the Fourier transform of the quasi-1D wavefunction as This relation is expected to give an accurate prediction of the dephasing time for all states where only 3 lattice sites are dominant.From this equation, one can see how the dephasing time tends to infinity in the limit of equally distributed chemical potentials, as observed in Fig. 2 of the main text.We can extend this to 5 sites, but it is not as trivial.One needs to numerically solve the tran- scendental equation for the smallest non-zero root t d .We compare the results from Eqs. (S12) and (S13) to the numerically obtained dephasing rate, γ = 1/t d , in Fig. S2(c), as presented in Fig. 2 of the main text, and find excellent agreement.

EXPERIMENTAL PROTOCOL
We prepare a 166 Er spin-polarized BEC similar to Ref. [4].The magnetic field during the evaporation is along the z-axis with absolute value |B| = B z = 1.9 G (a s = 80(1) a 0 ), see Fig. 1(a).The B-to-as conversion has been precisely mapped out in previous experiments [4,20].Before loading the lattice, we rotate the magnetic field direction along the y-axis in 50 ms and change its absolute value to set the scattering length.At this step, we typically achieve 5×10 4 atoms with more than 60 % condensed fraction in a cigar shape dipole trap with trapping frequencies ω x,y,z = 2π (240(3), 30(3), 217(1)) Hz.For our experiments, the atoms are then loaded in a 1D lattice by a 20 ms exponential ramp of the lattice depth.This is the experimental protocol used in Fig 1, 2, and 3(a).
To study the reversibility of the interaction-induced transition to a single lattice site (3(b)), i.e. the evolution of the contrast due to a change of the scattering length, we employ a different protocol from the one above.In fact, in our experiment, the magnetic field along the ydirection can be changed on a timescale of 20 ms, which is slower compared to the z-direction ( 1 ms).For this dataset, we prepare the BEC with B = (0, 0.25, 1) G and then we load the lattice as described above.We then linearly ramp the field in 20 ms to B = (0, 0.25, 0) G and record the time evolution.In Fig. 3(b), we study the contrast evolution after the ramp.For the black dataset, the magnetic field is quenched back to the initial value after 10 ms.
For Fig 4, we extract the atom number condensed in the lattice by releasing the cloud from the combined ODT-lattice trap and by performing an absorption imaging after 30 ms of TOF.We integrate the density along the lattice axis and use a double Gaussian fit on the integrated density profile.We repeat the sequence 4-8 times for every scattering length.At low scattering lengths, we find a decreased number of condensed atoms, see Fig. 4. We attribute this to an increase of three-body loss in the vicinity of a Feshbach resonance [4] and the increased density of the groundstate.min + 1. Below, (b) probability distribution of τ , given by numerical integration of A dAe −χ 2 (A,τ )/2 and normalization to 1. as=58.8 a0.The dashed line indicates the fit result and the shaded area the 68% confidence interval.

ANALYSIS OF MOMENTUM DISTRIBUTION DURING BLOCH OSCILLATION
When the Bloch oscillation dephases, the width of the momentum distribution increases with time [82].To evaluate the dephasing rate we analyze the 1D momentum distribution along z, n(q z ), as a function of the holding time.Because of our limited vertical optical access, the 1D lattice is not perfectly aligned with the z (gravity) direction.We measure a tilt of 9(1)°.Such a tilt effectively weakens the radial trapping strength, limiting our observation time to 12 ms, which anyhow allows us to observe up to 25 BO period.
From n(q z ), we can extract the maximum position (q max z ) and the quantity |q z | , given by This quantity is proportional to the width of the distribution.In Fig. S3, we report |q z | for a s = 65.7(1.0) a 0 .To quantify the dephasing rate γ, we apply a linear fit to |q z | .For the fit, we select only the points at the center of the Brillouin zone, up to the time when |q z | is reaching the fully dephased configuration, 0.5 k.Indeed, when the cloud is at the edge of the Brillouin zone, |q z | is artificially increased and it does not represent the dephasing, as shown in Fig. S3.We define the dephasing rate γ as the inverse of the time τ that the fitted function needs to reach the value 0.5 k.Thus, using the fit parametrization A(t−τ )+0.5, where A and τ are the fitting variables and t is the time, we can directly extract τ and its inverse γ.
To determine the uncertainties with our non-linear parametrization, we analyze the χ 2 (A, τ ).We estimate the uncertainties on our data points by assuming equal statistical fluctuations around our fitting model and using the expected value χ 2 = N data − 2. Figure S4 shows a clear asymmetric shape for χ 2 , indicating asymmetric uncertainties on our fit parameters.As we are only interested in the uncertainties on τ , we consider P χ 2 (τ ) = 1 N A dAe −χ 2 (A,τ )/2 , with N a normalization constant.P χ 2 (τ ) corresponds to the probability distribution of τ for our fitting model.Finally, from P χ 2 (τ ), we define the 68% confident interval of our dephasing rate γ shown in Fig. S4.
In order to compare our experimental data with the theoretical predictions, we repeat the same analysis with the data from the 1D discrete model.Since in the experiment the condensed atom number changes with the scattering length, see Fig. 4, the atom number consid-ered in the theoretical simulations varies accordingly.In Fig. 2, we account for the experimental fluctuations by taking an interval of ±20% of the BEC atoms number.For each scattering length, we determine the extreme values of γ in the ±20% range, which we use to create the shaded area.

CONTRAST OF THE INTERFERENCE PATTERN
The density modulation that usually characterizes a BEC loaded into a 1D lattice can be experimentally extracted from the matter-wave interferometry after a TOF expansion [59].To study the transition to one single occupied lattice site, we record the density distribution as a function of a s .In more details, for each picture we perform a Fourier transform (FT) of the integrated momentum distribution, n(q z ).In the contact dominated regime, the lattice induces two sidepeaks at ±q * z in n(q z ).Consequently, in the FT analysis the peaks are at z * λ lattice .The visibility of the interference pattern is then estimated as n FT (|z * |)/n FT (0).

FIG. 2 .
FIG.2.Dephasing rate and chemical potential distributions.(a) Experimental dephasing rate γ (circles) as a function of scattering length as.The green solid line shows the theory result, with an uncertainty region (shaded area) accounting for 20% atom number variation.The blue dashed line shows the theory expectation without LHY.The gray dotdashed line gives the prediction of the semi-analytic approximation for γ.Error bars show the 68% confidence interval[50].(b) Chemical potential per lattice site µj extracted from the discrete model for as = (59, 60, 65.5, 70) a0 (1, 2, 3, 4).The green area depicts the LHY contribution to µj.

FIG. 3 .
FIG. 3. Interaction-induced localization.(a) Contrast of the interference pattern after loading the lattice at different as.The green dot-dashed (black solid) line represents the result of the 1D discrete model (3D eGPE) multiplied by 0.7.The insets show the respective density distributions along z of the 1D discrete model (bars) and 3D eGPE (lines) and corresponding experimental averaged interference patterns after TOF expansion (1,2).(b) Dynamic evolution of the contrast quenching back (filled circles) or holding as (open circles); see text.The error bars represent the standard error on the mean over 4-6 repetitions.

FIG. 4 .
FIG. 4. Phase diagram and energy landscapes.(a) Phase diagram as a function of as and atom number.The white region denotes a trap-bound BEC extended over several lattice sites.The colored regions denote quasi-2D selfbound solutions: a droplet (green), a soliton (blue), each either extended over several lattice sites (lighter shade) or localized (darker shade, >95% of the atoms are localized in the central lattice plane).Circles show our experimental data points from Fig. 3(a).Inset (a), (b-c) Energy landscapes as a function of the radial widths lx and ly, in units of the radial harmonic oscillator lengths x ho = 0.50(1) µm and y ho = 1.42(1) µm, respectively, with (left) and without (right) the radial harmonic trap, for (inset (a)) BEC (as, N ) = (70 a0, 1.5 × 10 4 ), (b) droplet (as, N ) = (65 a0, 1.5 × 10 4 ) and (c) soliton (as, N ) = (51.5 a0, 0.4 × 10 4 ) regimes, with darker shading at the minima.(d) Radial width lx versus N for as = 51.5 a0.The dashed line indicates the soliton-to-droplet transition point, and the circles indicate the position of (b-c).
Supplemental materials: Strongly dipolar gases in a one-dimensional lattice: Bloch oscillations and matter-wave localization G. Natale, T. Bland, S. Gschwendtner, L. Lafforgue, D. S. Grün, A. Patscheider, M. J. Mark, and F. Ferlaino THEORETICAL MODEL FIG. S1.Interaction energy contributions.Scattering length dependency of the individual interaction contributions of the ground state solutions from the 1D model, calculated for N = 10 4 atoms.
S10)where w(k) is the momentum space Wannier function, and phases φ j are given above.If all interactions are set to zero this function is initially a delta function situated at k = 0 and moves in k-space as k = k − F ext t/ .Interactions broaden C(k, t), leading to a dephasing of coefficients c j .Fig.S2(a) depicts | C(k, t)| 2 as a function of k at different times t, normalized to | C(0, 0)| 2 .We extract an analytic approximation to the dephasing time by considering the temporal behaviour of the point | C(0, t)| 2 , i.e. at k = 0.During dephasing this point rapidly decreases through interference between neighboring sites.This quantity is plotted in Fig.S2(b) for a few example scattering lengths.It reaches the threshold α/ C(0, 0) = 0.5 at the dephasing time t = t d , where many k-modes are now highly occupied.This time can be found through the smallest positive solution of | C(0, t d )| 2 = α can be only found in limiting cases.For the three lattice site case, with j = −1, 0, 1 and noting the symmetry of |c j | = |c −j | we obtain

FIG. S2 .
FIG. S2.Analytic dephasing rate.(a) Evolution of the function C, with k normalized to the Brillouin zone in the moving frame, and as = 60.5 a0.Here, the solution of Eq. (S12) is t d = 0.59s.(b) Time evolution of the central point of C, showing when | C| 2 crosses α = 0.5.The function C is scaled to the value at C(0, 0).(c) Analytic dephasing rate (γ = 1/t d ) obtained for the 3 lattice site approximation Eq. (S12) and the 5 lattice site approximation Eq. (S13), compared to the numerically obtained value from a real-time simulation of the discrete model, Eq. (S4).
FIG. S3.Evolution of the |qz| .In the figure |qz| as a function of time for as=58.8a0.The points with red edges are the one selected for the fit.The red line corresponds to the fit result A(t − τ ) + 0.5

3 FIG
FIG. S4.Uncertainties χ 2 analysis.(a) χ 2 as a function of A and τ .The black dash lines correspond to the value χ 2 min + 1. Below, (b) probability distribution of τ , given by numerical integration of A dAe −χ 2 (A,τ )/2 and normalization to 1. as=58.8 a0.The dashed line indicates the fit result and the shaded area the 68% confidence interval.