Non-Equilibrium Mass Transport in the 1D Fermi-Hubbard Model
S. Scherg, T. Kohlert, J. Herbrych, J. Stolpp, P. Bordia, U. Schneider, F. Heidrich-Meisner, I. Bloch, M. Aidelsburger
NNon-Equilibrium Mass Transport in the 1D Fermi-Hubbard Model
S. Scherg,
1, 2
T. Kohlert,
1, 2
J. Herbrych,
3, 4
J. Stolpp,
1, 5, 6
P. Bordia,
1, 2
U. Schneider,
1, 2, 7
F. Heidrich-Meisner, I. Bloch,
1, 2 and M. Aidelsburger
1, 2 Fakult¨at f¨ur Physik, Ludwig-Maximilians-Universit¨at M¨unchen, Munich, Germany Max-Planck-Institut f¨ur Quantenoptik, 85748 Garching, Germany Department of Physics and Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USA Materials Science and Technology Division Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA Institute for Theoretical Physics, Georg-August-Universit¨at G¨ottingen, 37077 G¨ottingen, Germany Arnold Sommerfeld Center for Theoretical Physics,Ludwig-Maximilians-Universit¨at M¨unchen, 80333 Munich, Germany Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, UK
We experimentally and numerically investigate the sudden expansion of fermions in a homogeneousone-dimensional optical lattice. For initial states with an appreciable amount of doublons, we observea dynamical phase separation between rapidly expanding singlons and slow doublons remaining inthe trap center, realizing the key aspect of fermionic quantum distillation in the strongly-interactinglimit. For initial states without doublons, we find a reduced interaction dependence of the asymptoticexpansion speed compared to bosons, which is explained by the interaction energy produced in thequench.
Many-body physics in one dimension (1D) differs innumerous essential aspects from its higher-dimensionalcounterparts. Several familiar concepts, such as Fermi-liquid theory [1, 2], are not applicable in 1D. More-over, many 1D models are integrable, meaning thatthere exist exact solutions. Examples include the Lieb-Liniger model [3], the Heisenberg chain [4] or the 1DFermi-Hubbard model (FHM) [5]. These models ex-hibit extensive sets of conserved quantities that pre-vent thermalization [6–11] and can, in lattice systems,lead to anomalous transport properties [12–15]. Cold-atom experiments offer the possibility to study trans-port properties of strongly-correlated quantum gases ina clean environment. Their excellent controllability en-abled far-from-equilibrium experiments [16–20] as well asclose-to-equilibrium measurements in the linear-responseregime [21–24] both in extended lattices and mesoscopicsystems [25–27].Here, we investigate mass transport in the 1D FHMin far-from-equilibrium expansion experiments [18–20],where an initially trapped gas is suddenly released into ahomogeneous potential landscape as illustrated in Fig. 1.There are two distinct regimes of interest in sudden-expansion studies: the asymptotic one, where the ex-panding gas has become dilute and effectively non-interacting [28–36] and the transient regime, where thedynamical quasi-condensation of hardcore bosons [37–41]and quantum distillation [20, 42–44] have been found.Quantum distillation occurs for large interactions. Itrelies on the dynamical demixing of fast singlons (oneatom per site) and slow doublons (two atoms per site)during the expansion: while isolated doublons only movewith a small effective second-order tunneling matrix ele-ment J eff = 2 J /U (cid:28) J for U (cid:29) J [46, 47], neighboringsinglons and doublons can exchange their positions viafast, resonant first-order tunneling processes. Thus, af-ter opening the trap, singlons escape from regions of thecloud initially occupied by singlons and doublons, leading (a) Initial state with doublons (b)
Initial state without doublons x dU J x JJ
10 20 30 40 50Site index i T i m e t ( ¿ )
10 20 30 40 50Site index i ̂ n i ® FIG. 1.
Schematics of the expansion experiment.
Top: Initial state of the harmonically trapped two-componentFermi gas with (a) singlons (red) and doublons (blue) and (b)only singlons in an optical lattice. After quenching to lowerlattice depths and removing the harmonic trap, fermions ex-pand in a homogeneous 1D lattice; J = h · . U refers to the on-site interactionstrength and d is the lattice constant. The expansion dynam-ics is dominated by first-order processes: (a) the resonantexchange of singlon and doublon positions leads to quantumdistillation, (b) the dynamical formation of doublons resultsin reduced asymptotic expansion velocities. Bottom: time-dependent density-matrix renormalization group (tDMRG)simulations of the atomic density (cid:104) ˆ n i (cid:105) for U = 20 J as a func-tion of time t in units of the tunneling time τ = (cid:126) /J . to a spatial separation of the two components. Withoutan increase of the doublon density in the central region,this regime is termed weak quantum distillation [44].Ideal initial-state conditions can lead to a strong versionof quantum distillation, where the spatial separation of a r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p Expansion time t ( ¿ ) R s ( d ) , R d ( d ) U =5 J singlonsdoublons 0 10 20 30 40 Expansion time t ( ¿ ) U =20 J
150 0 1500.00.10.20.3 ½ s , ½ d ( a . u . )
150 0 1500 10 20 30 40
Expansion time t ( ¿ ) N c d / N c s t ( ¿ ) N s , N d t ( ¿ ) N s , N d (a) (b) (c)(d) U =20 Jt =0 t =40 ¿ Position x ( d )
40 d x
FIG. 2.
Dynamical phase separation of singlons and doublons.
Half-width-at-half-maximum (HWHM) size R s,d ofsinglon (s, red) and doublon (d, blue) clouds as a function of time for (a) U = 5 J and (b) U = 20 J . The dashed lines illustratethe hypothetical expansion of a non-interacting doublon cloud with effective tunneling J eff [45]. Insets: Number of atoms onsingly- and doubly-occupied sites, N s and N d , as a function of time. (c) Experimental snapshots of the integrated line densitiesfor singlon and doublon clouds, ρ s ( x ) and ρ d ( x ), at t = 0 (left) and t = 40 τ (right). (d) Ratio N cd /N cs of atom numbers ondoubly- and singly-occupied sites in the central region of the cloud (red rectangle in the inset) as a function of time for U = 20 J .Every data point is averaged over four measurements and error bars represent the standard-error-of-the-mean. Solid lines areguides to the eye. singlons and doublons yields a contraction of the dou-blon cloud radius [Fig. 1(a)]. In the extreme limit thisdistillation could be used to purify a finite-temperatureband insulator [42], thereby dynamically generating low-entropy regions. This would represent a major advance-ment in the ongoing quest of realizing fermionic many-body physics at the lowest entropy scales [48–52]. So far,experimental evidence for weak quantum distillation hasonly been found for bosons [20] at intermediate interac-tion strengths, where doublons can decay into singlonson time scales relevant to quantum distillation.In this work, we investigate the non-equilibrium masstransport in the 1D FHM starting from initial productstates in deep optical lattices [18, 19], while close-to-thermal initial states were used in Ref. [20]. The expan-sion dynamics is initiated by two simultaneous quenches:a sudden increase of the tunnel coupling, resulting in aquench from almost infinite to finite U/J , and a sud-den removal of the harmonic trap (Fig. 1). We prepareinitial product states with or without doublons (Fig. 1)and quantitatively investigate the time evolution of sin-glon and doublon densities individually. For initial stateswith doublons, we find a distinct dynamical phase sepa-ration between singlons and doublons, which is the fun-damental mechanism of fermionic quantum distillation[Fig. 1(a)]. For initial states without doublons [Fig. 1(b)],we study interaction effects in the asymptotic expansionvelocities [35]. We observe that the cloud expands rapidlyat all interaction strengths with slightly smaller velocitiesat intermediate values, in agreement with our numericalsimulations. This can be interpreted in terms of the in-teraction energy produced in the quench of
U/J , which leads to the dynamical formation of doublons [19, 33, 53].
Experiment.
We prepare a degenerate Fermi gas of30(1) × K atoms in a crossed dipole trap at the ini-tial temperature
T /T F = 0 . T F is the Fermitemperature. The gas consists of an equal mixture oftwo spin components corresponding to the states |↑(cid:105) = | m F = − / (cid:105) and |↓(cid:105) = | m F = − / (cid:105) in the F = 9 / λ x = 532 nm and latticeconstant d = λ x / x direction and λ ⊥ = 738 nmin the transverse directions. While the main lattice along x is initially loaded to 20 E rx , the transverse lattices aresimultaneously ramped to a depth of 33 E r ⊥ , where theyremain during the whole sequence to realize individual1D systems. Here, E rj = (cid:126) k j / (2 m ) are the respectiverecoil energies with j ∈ { x, ⊥} , k j = 2 π/λ j denotes thecorresponding wave vector and m is the mass of K.Holding the atoms in the deep initial lattice for 20 ms de-phases remaining correlations between neighboring sites,such that the resulting state can be approximated as aproduct state | ψ (cid:105) = (cid:81) i ∈ trap (cid:16) ˆ c † i ↑ (cid:17) n i ↑ (cid:16) ˆ c † i ↓ (cid:17) n i ↓ | (cid:105) , whereˆ c † iσ is the fermionic creation operator, n iσ ∈ { , } , σ ∈ {↑ , ↓} and i is the lattice-site index. The spin orien-tations are expected to be distributed randomly amongthe sites and the average number of atoms per lattice sitein the center of the cloud (cid:104) ˆ n i (cid:105) = (cid:80) σ (cid:104) ˆ n iσ (cid:105) is estimatedto be (cid:104) ˆ n i (cid:105) (cid:46) . n iσ = ˆ c † iσ ˆ c iσ is the density opera-tor. The fraction of atoms on doubly-occupied sites n d = N d / ( N s + N d ) in the initial state can be tuned via the in-teraction strength during the loading process employinga Feshbach resonance at 202 . N s ( N d ) de-notes the number of particles on singly(doubly)-occupiedsites. The dynamics starts with suddenly quenching themain lattice to 8 E rx . Simultaneously, the strength of thedipole trap is adjusted to compensate the anti-confiningharmonic potential introduced by the lattice [45]. Oursystem is then described by the homogeneous 1D FHMˆ H = − J (cid:88) i,σ = ↑ , ↓ (cid:16) ˆ c † iσ ˆ c i +1 σ + h.c. (cid:17) + U (cid:88) i ˆ n i ↑ ˆ n i ↓ . (1)After a variable expansion time t the on-site population isfrozen by suddenly increasing the lattice depth to 20 E rx .Subsequently, the cloud is imaged in-situ using high-fieldimaging either with or without removing doublons [45].By combining these images, the dynamics of singlons anddoublons can be resolved individually. Quantum distillation.
We characterize the dynam-ics by monitoring the singlon and doublon clouds as afunction of the expansion time for an initial state with n d = 0 . U (cid:29) W , W = 4 J [54], since in this case the interaction energyreleased in the doublon decay cannot be transferred intokinetic energy of singlons in low-order processes [46, 54].This is in agreement with our observations [insets inFigs. 2(a), (b)], where for U = 5 J we witness a fast dou-blon decay of about 25% in the early stages of the ex-pansion t (cid:46) τ , which is accompanied by a compatibleincrease of the singlon number. In contrast, both num-bers remain approximately constant for U = 20 J . Exceptfor a small residual decay, which is attributed to light-assisted losses of doublons [55], this enables us to probethe dynamical phase separation of singlons and doublonsat approximately constant doublon numbers.We study the phase separation by extracting the cloudsizes R s,d ( t ) at half-width-at-half maximum (HWHM).We observe a rapidly-expanding singlon cloud, which hasapproximately doubled in size at t = 40 τ . In contrast, thedoublon cloud size grows much slower and we even ob-serve a weak shrinking of the cloud for U = 20 J . Forcomparison, we show the expected expansion of a fic-titious cloud of non-interacting doublons expanding ac-cording to J eff [47]. The difference highlights the non-trivial nature of this transient dynamics. The dynamicalphase separation is even more evident in the comparisonof the integrated line densities of singlons and doublons at t = 0 and t = 40 τ for our strongest interactions [Fig. 2(c)].Clearly, the singlons expand significantly, while the dou-blons essentially remain in the center of the cloud. As aconsequence, the ratio of atom numbers on doubly- andsingly-occupied sites N cd /N cs in the center of the cloud in-creases by about 40% [Fig. 2(d)]. While this signal couldin principle be caused by 1D systems with a low doublonfraction, a quantitative analysis based on our measuredinitial density distributions shows that their contributionto the signal is negligible [45]. Hence, our data establish Expansion time t ( ¿ ) −0.6−0.5−0.4−0.3−0.2−0.10.00.1 ¢ R d n =1.25 n =1.00 n =0.83 5 10 15 20 U / J −0.10.00.10.2 ¢ R d t =40 ¿ FIG. 3. tDMRG results for the relative doublon cloudsize.
Main panel: Relative doublon could size ∆ R d ( t ) forthree different initial uniform densities n at U/J = 20. Theinitial state consists of 12 singlons, 4 doublons and { , , } holons, respectively [45]. The solid lines end at time t max ,when the width of the singlon cloud increased to ∆ R s = 0 . t = 40 τ .Inset: Experimental data for ∆ R d as a function of the inter-action strength at t = 40 τ , which was evaluated using linearfits to the time traces R d ( t ) as shown in Figs. 2(a), (b) [45]. clear evidence for fermionic quantum distillation in theweak regime in a non-equilibrium mass-transport exper-iment. tDMRG results for transient dynamics. Quantum dis-tillation in the strong regime can further lead to a shrink-ing of the doublon cloud. The precise amount depends onthe number of singlons initially confined in the doubloncloud, the initial density and the cloud size [42, 44, 45].Here, we focus on the role of the initial density, whichhas the largest influence. Figure 3 shows tDMRG sim-ulations of the relative change of the doublon cloud size∆ R d ( t ) = R d ( t ) /R d (0) − n = ( N s + N d ) /L init in anideal box trap of length L init for constant n d = 0 . R d indicate a shrinking of the dou-blon cloud, while ∆ R d > n = 1 .
25, we observe a large decrease of ∆ R d ( t ).This effect is substantially reduced for smaller densities(Fig. 3) due to the presence of holons (empty sites), whichremain trapped between doublons on the time scales ofthe quantum distillation process [44]. Additionally, thedynamics becomes slower, both due to holons and due tothe larger cloud sizes used for simulations with smalleraverage densities [44]. Despite these quantitative differ-ences, however, we find that the fundamental aspect ofquantum distillation, i.e., the dynamical phase separa-tion of singlons and doublons, is generally robust.For comparison, we show the experimentally measuredrelative changes ∆ R d as a function of interaction strength U / J v r ( d / ¿ ) ExperimenttDMRG fermionstDMRG N ́eeltDMRG bosons
FIG. 4.
Radial expansion velocities v r . Experiment(circles) and tDMRG simulations for fermions (red dark-shaded diamonds) and bosons (green light-shaded diamonds,from [19]) as a function of
U/J . Solid lines are guides tothe eye. The grey dashed line indicates v r = √ d/τ in thelimiting cases U/J = 0 and
U/J → ∞ . All initial states inthe numerical simulations have a uniform average density of n = 1 in a box with initial size L init = 10. [inset in Fig. 3]. For all interactions the time traces of thedoublon HWHM are fitted with a linear function to cal-culate ∆ R d at the maximum expansion time t = 40 τ [45].We observe that ∆ R d (40 τ ) approaches zero with increas-ing interaction strength and becomes slightly negative at U/J = 20. In order to facilitate a comparison betweenexperiment and numerics, where much smaller particlenumbers are used, we define a time t max for the simu-lation at which the relative singlon cloud size ∆ R s hasreached the same value as in the experiment (Fig. 3).Our numerical results indicate that the contraction isnot completed at this time. In the experiment this timeis limited by the degree of flatness of the homogeneouspotential. The remaining difference between the numer-ical and experimental results is most likely due to otherinitial-state properties, such as inhomogeneous densitydistributions and the averaging over several 1D systemswith different initial-state properties [45]. Asymptotic dynamics and interaction effects.
Here, wefocus on the dynamics of the whole cloud for initial stateswith a negligible doublon fraction ( n d < . r = (cid:80) i (cid:104) ˆ n i (cid:105) ( i − i ) d / ( N s + N d )of the time-dependent density distribution (see [45] fordetails on the analysis), which is routinely computed innumerical simulations [32, 33, 35]; here i is the center-of-mass of the initially trapped gas. From the time de-pendence of r , we extract the asymptotic radial velocity v r by fitting √ r = (cid:112) r + v r t , where r is the initialsize of the cloud [45]. Figure 4 shows v r as a function of U/J . We find v r = 1 . d/τ for U = 0 and U = 20 J ,whereas for intermediate interactions U ∼ J , the ra-dial velocity decreases weakly. Note that for U (cid:29) W , the mass transport in the 1D FHM in the absence ofdoublons becomes identical to a non-interacting gas ofspinless fermions and thus it behaves exactly like hard-core bosons in 1D with v r = √ d/τ [19]. The valuesin the limiting cases agree with these theoretical predic-tions for free fermions expanding from our initial state.Remarkably, compared to the Bose-Hubbard model [19],the reduction of v r at intermediate interaction strengthsis much weaker (Fig. 4).Starting from the limit of very strong interactions, theinteraction dependence of v r can be understood in a two-component picture of independent singlon and doublongases [56, 57]: The dynamically generated doublons un-dergo a quantum distillation mechanism and are theninert on the time scales of the experiment. Thus, themore doublons are generated, the less kinetic energy isavailable for the rapidly expanding singlons. Focusing onthe quantitative difference between the v r ( U ) curves forbosons and fermions, which is the main result of the datapresented in Fig. 4, two aspects are important. First,in the case of fermions, doublons can only be gener-ated between sites with fermions of different spin orienta-tion [33]. The initial state that has the most of such ↑ - ↓ neighbors is the N´eel state, and this initial state leadsto the most pronounced minimum of v r (Fig. 4, [33]).In order to compare to the experiment, we average overmany 1D systems with random spin orientations for abalanced spin mixture (dark red diamonds in Fig. 4).This averaging leads to a weaker minimum in v r than forthe N´eel state and is in agreement with our experimentaldata. The second reason for the stronger minimum in v r for bosons is the fact that the interaction energy canbecome much larger, since larger local occupancies arepossible [19]. In order to test whether the observed v r can primarily be understood as a function of the inter-action energy in the system after the formation of dou-blons, we show data for different U/J and different spinconfigurations versus interaction energy in the Supple-mental Material (Fig. S8 in [45]). The data for bosonslie well outside the accessible range of interaction ener-gies for fermions because of higher site occupations, butfall onto an extrapolation of the fermionic data. Hence,the integrability of the 1D FHM does not seem to be thedominant reason for the differences to the bosonic case.An interesting extension would be the calculation of ex-pansion velocities by exploiting the integrability alongthe lines of [35, 58], which we leave for future work.
Summary and Outlook.
We investigated the suddenexpansion of an interacting cloud of fermions. Startingfrom an initial product state with an appreciable dou-blon fraction, we observed a dynamical phase separationbetween singlons and doublons, theoretically known asfermionic quantum distillation in the weak regime. Ad-ditionally, we analyzed radial velocities for different in-teraction strengths using initial states consisting purelyof singlons. We found a decrease of the radial velocitiesat weak interactions and attributed this effect to dynam-ically generated doublons. The weak decrease of radialvelocities of expanding fermions compared to bosons isdue to the Pauli principle leading to a crucial depen-dence of the radial velocities on the initial spin config-uration. Future experiments could use the singlon anddoublon resolved scheme to detect signatures of FFLOstates [58–60] in the expansion velocity of the unpairedspin component. Moreover, it would be intriguing to ob-serve the strong version of quantum distillation, resultingin the dynamical formation of low-entropy regions. Thiscould be achieved by optimizing the initial-state proper-ties and improved imaging techniques, such as microwavedressing to isolate central 1D systems, where the condi- tions for strong quantum distillation are best, or usingquantum-gas microscopes [61, 62].
Acknowledgments.
We acknowledge financial sup-port by the European Commission (UQUAM grant no.319278, AQuS) and the Nanosystems Initiative Munich(NIM) grant no. EXC4. J.S. and F.H.-M. were sup-ported by the DFG (Deutsche Forschungsgemeinschaft)Research Unit FOR 1807 under grant no. HE 5242/3-2and by startup funds via SFB 1073 at the University ofG¨ottingen. J.H. was supported by the US Departmentof Energy (DOE), Office of Basic Energy Sciences (BES),Materials Sciences and Engineering Division. [1] G. Baym and C. Pethick, “Landau fermi-liquid theory,”Wiley, New York (1991).[2] T. Giamarchi, “Quantum physics in one dimension,” Ox-ford university press (2004).[3] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac,and M. Rigol, “One dimensional bosons: From condensedmatter systems to ultracold gases,” Rev. Mod. Phys. ,1405 (2011).[4] A. Kl¨umper, “Integrability of quantum spin chains: The-ory and applications to the spin-1/2 XXZ chain,” Lect.Notes Phys. , 349 (2004).[5] F.H.L. Essler, H. Frahm, F. G¨ohmann, A. Kl¨umper, andV. E. Korepin, “The one-dimensional Hubbard model,”Cambridge University Press (2005).[6] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, “Re-laxation in a completely integrable many-body quantumsystem: An ab initio study of the dynamics of the highlyexcited states of 1d lattice hard-core bosons,” Phys. Rev.Lett. , 050405 (2007).[7] L. Vidmar and M. Rigol, “Generalized Gibbs ensemblein integrable lattice models,” J. Stat. Mech. Theor. Exp. , 064007 (2016).[8] T. Kinoshita, T. Wenger, and D. S. Weiss, “A quantumNewton’s cradle,” Nature , 900 (2006).[9] S. Hofferberth, I. Lesanovsky, B. Fisher, T. Schumm,and J. Schmiedmayer, “Non-equilibrium coherence dy-namics in one-dimensional Bose gases,” Nature , 324(2007).[10] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa,B. Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler,and J. Schmiedmayer, “Relaxation and prethermaliza-tion in an isolated quantum system,” Science , 1318–1322 (2012).[11] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler,M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer,and J. Schmiedmayer, “Experimental observation of ageneralized Gibbs ensemble,” Science , 207 (2015).[12] X. Zotos, F. Naef, and P. Prelovˇsek, “Transport andconservation laws,” Phys. Rev. B , 11029 (1997).[13] F. Heidrich-Meisner, A. Honecker, and W. Brenig,“Transport in quasi one-dimensional spin-1/2 systems,”Eur. Phys. J. Special Topics , 135 (2007).[14] R. Vasseur and J.E. Moore, “Nonequilibrium quantumdynamics and transport: from integrability to many-bodylocalization,” J. Stat. Mech. Theor. Exp. , 064010(2016). [15] C. Karrasch, T. Prosen, and F. Heidrich-Meisner, “Pro-posal for measuring the finite-temperature drude weightof integrable systems,” Phys. Rev. B , 060406 (2017).[16] H. Ott, E. de Mirandes, F. Ferlaino, G. Roati, G. Mod-ugno, and M. Inguscio, “Collisionally induced trans-port in periodic potentials,” Phys. Rev. Lett. , 160601(2004).[17] N. Strohmaier, Y. Takasu, K. G¨unter, R. J¨ordens,M. K¨ohl, H. Moritz, and T. Esslinger, “Interaction-controlled transport of an ultracold fermi gas,” Phys.Rev. Lett. , 220601 (2007).[18] U. Schneider, L. Hackerm¨uller, J. P. Ronzheimer, S. Will,S. Braun, T. Best, I. Bloch, E. Demler, S. Mandt,D. Rasch, and A. Rosch, “Fermionic transport andout-of-equilibrium dynamics in a homogeneous Hubbardmodel with ultracold atoms,” Nat. Phys. , 213 (2012).[19] J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodg-man, S. Langer, I. P. McCulloch, F. Heidrich-Meisner,I. Bloch, and U. Schneider, “Expansion dynamics of in-teracting bosons in homogeneous lattices in one and twodimensions,” Phys. Rev. Lett. , 205301 (2013).[20] L. Xia, J. Carrasquilla, A. Reinhard, J. M. Wilson,M. Rigol, and D. S. Weiss, “Quantum distillation andconfinement of vacancies in a doublon sea,” Nat. Phys. , 316–320 (2015).[21] W. Xu, W. Mcgehee, W. Morong, and B. DeMarco, “Badmetal in a fermi lattice gas,” arXiv:1606.06669 (2016).[22] R. Anderson, F. Wang, P. Xu, V. Venu, S. Trotzky,F. Chevy, and J. H. Thywissen, “Optical conductivityof a quantum gas,” arXiv:1712.09965 (2017).[23] P. T. Brown, D. Mitra, E. Guardado-Sanchez,R. Nourafkan, A. Reymbaut, S. Bergeron, A.-M. S. Trem-blay, J. Kokalj, D.A. Huse, P. Schauss, and W. S. Bakr,“Bad metallic transport in a cold atom Fermi-Hubbardsystem,” arXiv:1802.09456 (2018).[24] M. A. Nichols, L.W. Cheuk, M. Okan, T. R. Hartke,E. Mendez, T. Senthil, E. Khatami, H. Zhang, and M.W.Zwierlein, “Spin transport in a Mott insulator of ultra-cold fermions,” arXiv:1802.10018 (2018).[25] J.-P. Brantut, J. Meineke, D. Stadler, S. Krinner, andT. Esslinger, “Conduction of ultracold fermions througha mesoscopic channel,” Science , 1069–1071 (2012).[26] G. Valtolina, A. Burchianti, A. Amico, E. Neri, K. Xhani,J. A. Seman, A. Trombettoni, A. Smerzi, M. Zaccanti,M. Inguscio, and G. Roati, “Josephson effect in fermionicsuperfluids across the bec-bcs crossover,” Science , , 011053 (2018).[28] M. Rigol and A. Muramatsu, “Fermionization in an ex-panding 1d gas of hard-core bosons,” Phys. Rev. Lett. , 240403 (2005).[29] A. Minguzzi and D. M. Gangardt, “Exact coherent statesof a harmonically confined Tonks-Girardeau gas,” Phys.Rev. Lett. , 240404 (2005).[30] A. del Campo and J.G. Muga, “Dynamics of a Tonks-Girardeau gas released from a hard-wall trap,” EPL ,965 (2006).[31] D. Juki´c, B. Klajn, and H. Buljan, “Momentum dis-tribution of a freely expanding Lieb-Liniger gas,” Phys.Rev. A , 033612 (2009).[32] S. Langer, M. J. A. Schuetz, I. P. McCulloch,U. Schollw¨ock, and F. Heidrich-Meisner, “Expansion ve-locity of a one-dimensional, two-component Fermi gasduring the sudden expansion in the ballistic regime,”Phys. Rev. A , 043618 (2012).[33] L. Vidmar, S. Langer, I. P. McCulloch, U. Schneider,U. Schollw¨ock, and F. Heidrich-Meisner, “Sudden ex-pansion of Mott insulators in one dimension,” Phys. Rev.B , 235117 (2013).[34] A. S. Campbell, D. M. Gangardt, and K. V.Kheruntsyan, “Sudden expansion of a one-dimensionalBose gas from power-law traps,” Phys. Rev. Lett. ,125302 (2015).[35] Z. Mei, L. Vidmar, F. Heidrich-Meisner, and C. J.Bolech, “Unveiling hidden structure of many-body wavefunctions of integrable systems via sudden-expansion ex-periments,” Phys. Rev. A , 021607 (2016).[36] N. Schl¨unzen, S. Hermanns, M. Bonitz, and C. Verdozzi,“Dynamics of strongly correlated fermions: Ab initio re-sults for two and three dimensions,” Phys. Rev. B ,035107 (2016).[37] M. Rigol and A. Muramatsu, “Emergence of quasicon-densates of hard-core bosons at finite momentum,” Phys.Rev. Lett. , 230404 (2004).[38] I. Hen and M. Rigol, “Strongly interacting atom lasersin three-dimensional optical lattices,” Phys. Rev. Lett. , 180401 (2010).[39] M. Jreissaty, J. Carrasquilla, F.A. Wolf, and M. Rigol,“Expansion of Bose-Hubbard Mott insulators in opticallattices,” Phys. Rev. A , 043610 (2011).[40] L. Vidmar, J. P. Ronzheimer, M. Schreiber, S. Braun,S. S. Hodgman, S. Langer, F. Heidrich-Meisner, I. Bloch,and U. Schneider, “Dynamical quasicondensation ofhard-core bosons at finite momenta,” Phys. Rev. Lett. , 175301 (2015).[41] L. Vidmar, D. Iyer, and M. Rigol, “Emergent eigen-state solution to quantum dynamics far from equilib-rium,” Phys. Rev. X , 021012 (2017).[42] F. Heidrich-Meisner, S. R. Manmana, M. Rigol, A. Mu-ramatsu, A. E. Feiguin, and E. Dagotto, “Quantum dis-tillation: Dynamical generation of low-entropy states ofstrongly correlated fermions in an optical lattice,” Phys.Rev. A , 041603 (2009).[43] D. Muth, D. Petrosyan, and M. Fleischhauer, “Dynamicsand evaporation of defects in Mott-insulating clusters ofboson pairs,” Phys. Rev. A , 013615 (2012).[44] J. Herbrych, A. E. Feiguin, E. Dagotto, and F. Heidrich- Meisner, “Efficiency of fermionic quantum distillation,”Phys. Rev. A , 033617 (2017).[45] See the Supplemental Material, which includes Refs. [63–71], not cited in the main text, for details on the initialstate preparation, the creation of a flat potential, the con-trol of the doublon fraction, the singlon-doublon resolvedread-out, the estimation of tubes with few doublons, thefits to the time traces of R d , as well as the extraction ofthe radial velocities. Additionally, details of the tDMRGsimulations are explained.[46] K. Winkler, G. Thalhammer, F. Lang, R. Grimm,J. Hecker Denschlag, A. J. Daley, A. Kantian, H. P.Bahler, and P. Zoller, “Repulsively bound atom pairsin an optical lattice,” Nature , 853 (2006).[47] R. Rausch and M. Potthoff, “Filling-dependent doublondynamics in the one-dimensional Hubbard model,” Phys.Rev. B , 045152 (2017).[48] J.-S. Bernier, C. Kollath, A. Georges, L. De Leo, F. Ger-bier, C. Salomon, and M. K¨ohl, “Cooling fermionicatoms in optical lattices by shaping the confinement,”Phys. Rev. A , 061601 (2009).[49] T.-L. Ho and Q. Zhou, “Universal cooling scheme forquantum simulation,” arXiv:0911.5506 (2009).[50] A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons,M. Kanasz-Nagy, R. Schmidt, F. Grusdt, E. Demler,D. Greif, and M. Greiner, “A cold-atom Fermi-Hubbardantiferromagnet,” Nature , 462 (2017).[51] C. S. Chiu, G. Ji, A. Mazurenko, D. Greif, andM. Greiner, “Quantum state engineering of a Hubbardsystem with ultracold fermions,” Phys. Rev. Lett. ,243201 (2018).[52] D. McKay and B. DeMarco, “Cooling in strongly cor-related optical lattices: prospects and challenges,” Rep.Prog. Phys. , 054401 (2011).[53] F. Heidrich-Meisner, M. Rigol, A. Muramatsu, A. E.Feiguin, and E. Dagotto, “Ground-state reference sys-tems for expanding correlated fermions in one dimen-sion,” Phys. Rev. A , 013620 (2008).[54] A. Rosch, D. Rasch, B. Binz, and M. Vojta, “Metastablesuperfluidity of repulsive fermionic atoms in optical lat-tices,” Phys. Rev. Lett. , 265301 (2008).[55] M. T. DePue, C. McCormick, S. L. Winoto, S. Oliver,and D. S. Weiss, “Unity occupation of sites in a 3d opticallattice,” Phys. Rev. Lett. , 2262 (1999).[56] J. Kajala, F. Massel, and P. T¨orm¨a, “Expansion dynam-ics in the one-dimensional Fermi-Hubbard model,” Phys.Rev. Lett. , 206401 (2011).[57] S. Sorg, L. Vidmar, L. Pollet, and F. Heidrich-Meisner,“Relaxation and thermalization in the one-dimensionalBose-Hubbard model: A case study for the interactionquantum quench from the atomic limit,” Phys. Rev. A , 033606 (2014).[58] C. J. Bolech, F. Heidrich-Meisner, S. Langer, I. P. Mc-Culloch, G. Orso, and M. Rigol, “Long-time behavior ofthe momentum distribution during the sudden expansionof a spin-imbalanced Fermi gas in one dimension,” Phys.Rev. Lett. , 110602 (2012).[59] J. Kajala, F. Massel, and P. T¨orm¨a, “Expansion dynam-ics of the Fulde-Ferrell-Larkin-Ovchinnikov state,” Phys.Rev. A , 041601 (2011).[60] H. Lu, L. O. Baksmaty, C. J. Bolech, and H. Pu, “Ex-pansion of 1d polarized superfluids: The Fulde-Ferrell-Larkin-Ovchinnikov state reveals itself,” Phys. Rev. Lett. , 225302 (2012). [61] E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaude-cerf, G. D. Bruce, and S. Kuhr, “Quantum-gas mi-croscope for fermionic atoms,” Nature Physics , 738(2015).[62] L. W. Cheuk, M. A. Nichols, M. Okan, T. Gersdorf, V. V.Ramasesh, W. S. Bakr, T. Lompe, and M. W. Zwierlein,“Single-atom imaging of fermions in a quantum-gas mi-croscope,” Phys. Rev. Lett. , 193001 (2015).[63] U. Schneider, L. Hackerm¨uller, S. Will, Th. Best,I. Bloch, T. A. Costi, R. W. Helmes, D. Rasch, andA. Rosch, “Metallic and insulating phases of repulsivelyinteracting fermions in a 3d optical lattice,” Science ,1520 (2008).[64] S. R. White and A. E. Feiguin, “Real-time evolutionusing the density matrix renormalization group,” Phys.Rev. Lett. , 076401 (2004).[65] A.J. Daley, C. Kollath, U. Schollw¨ock, and G. Vidal,“Time-dependent density-matrix renormalization-groupusing adaptive effective Hilbert spaces,” J. Stat. Mech.Theor. Exp. , P04005 (2004).[66] G. Vidal, “Efficient simulation of one-dimensional quan-tum many-body systems,” Phys. Rev. Lett. , 040502(2004).[67] U. Schollw¨ock, “The density-matrix renormalizationgroup in the age of matrix product states,” Annals ofPhysics , 96 (2011).[68] C. Boschi, Degli E., E. Ercolessi, L. Ferrari, P. Naldesi,F. Ortolani, and L. Taddia, “Bound states and expan-sion dynamics of interacting bosons on a one-dimensionallattice,” Phys. Rev. A , 043606 (2014).[69] T. Enss and J. Sirker, “Light cone renormalization andquantum quenches in one-dimensional Hubbard models,”New J. Phys. , 023008 (2012).[70] A. Bauer, F. Dorfner, and F. Heidrich-Meisner, “Tem-poral decay of N´eel order in the one-dimensional Fermi-Hubbard model,” Phys. Rev. A , 053628 (2015).[71] F. Dorfner, “Doctoral thesis,” Ludwig-Maximilians Uni-versity Munich (2016). SUPPLEMENTAL MATERIAL
S1. Initial state preparation
1. Experimental sequence
The initial state preparation for the measurements inthe main paper starts with loading the atoms into thethree-dimensional (3D) optical lattice using a sequenceof different linear ramps. First, within 7 ms, the lat-tice along the x direction and the transversal latticesare ramped to a depth of 1 E rx and 1 E r ⊥ , respectively.Then, after waiting for 100 ms, the depth of all three lat-tices is further increased to 8 E rx and 8 E r ⊥ , respectively,during 75 ms. Finally, the lattice along the x directionis ramped up to 20 E rx and the transversal lattices reachtheir final depth of 33 E r ⊥ within 15 ms, freezing the pop-ulations of singlons and doublons. Afterwards, an addi-tional superlattice along the x direction is added at adepth of 20 E s , where E s = h / (2 mλ s ) is the recoil en-ergy of the superlattice with wavelength λ s = 1064 nm.The phase of the superlattice was set so as to create tilteddouble wells along the x direction, in order to decreaseresidual dynamics and remaining correlations betweenneighboring sites.While holding the atoms in the deep 3D optical latticefor 25 ms, both the dipole trap strength and the mag-netic field strength are adjusted to their target valuesduring the expansion of the cloud. We ramp the mag-netic field within 15 ms to change the scattering lengthfrom a s = − a (attractive loading to generate initialstates with doublons) or a s = 140 a (repulsive loading torealize initial states without doublons) to the scatteringlength which sets the desired interaction strength duringthe expansion of the cloud in the lattice (see Sec. S2).Moreover, the dipole traps along the x and the y direc-tion [trap frequency of ω x = ω y = ω = 2 π × E rx , E r ⊥ , E r ⊥ ) deep lattice dur-ing loading] are ramped to zero within 22 ms, whereasthe dipole trap along the z direction [trap frequency of ω z = 2 π × E rx , E r ⊥ , E r ⊥ )deep lattice during loading] is ramped within 22 ms suchthat the lattice potential is flat during the expansion (seeSec. S3). The preparation sequence terminates by re-moving the superlattice and quenching the lattice alongthe x direction within 10 µ s from 20 E rx to 8 E rx , in thisway initiating the expansion of the cloud.
2. Characterization
We characterize the initial state by estimating therenormalized cloud size R sc = R/ ( γ y γ z N σ ) / and the di-mensionless compression E t / (12 J ) during loading, where E t = V t [3 γ y γ z N σ / (4 π )] / and V t = mω d /
2. Here, γ y = 738 /
532 takes into account the different lattice con-stants between the x direction and the transversal y and −0.20.00.20.40.6 n d −250 −200 −150 −100 −50 0 50 100 a s ( a ) C l o u d s i z e ( d ) (a)(b) FIG. S1.
Control of the doublon fraction: (a) Fraction ofatoms on doubly-occupied sites n d depending on the loadingscattering length in units of the Bohr radius a and (b) theresulting cloud size in units of the lattice constant d . Toprepare initial states with a large number of doublons, we usea loading scattering length of a s = − a (dashed verticalline), which yields a significant fraction of atoms on doubly-occupied sites n d = 0 . z directions, γ z = ω z /ω = 184 /
54 takes into accountthe different harmonic confinement along the z directioncompared to x and y directions, N σ = 15(1) × isthe number of atoms per spin state σ , m is the massof K and d is the lattice constant along the x direc-tion. The parameters R sc and E t were previously used todistinguish between the Mott-insulating and the metal-lic regime of the Fermi-Hubbard model for repulsively-interacting fermions [63]. The renormalized cloud sizeis a rough estimate for the distance between two par-ticles in the same spin state and only depends on thedimensionless compression, the interaction strength andthe entropy set by the temperature in the pure har-monic trap. The dimensionless compression can be un-derstood as the ratio between the characteristic trap en-ergy E t , which is the Fermi energy of a non-interactinggas in the zero-tunneling limit, and the bandwidth 12 J in a 3D optical lattice. We estimate R sc = 0 . d and E t / (12 J ) = 0 . E rx , E r ⊥ , E r ⊥ ) deep lat-tice. This means that we work in the metallic regimewith a mean density in the center of the cloud of (cid:104) ˆ n i (cid:105) < S2. Control of the doublon fraction
The doublon fraction of the initial state can be set byproperly adjusting the scattering length a s during theloading of the atoms into the 3D optical lattice. While2large repulsive interactions ( a s >
0) result in initial stateswith a negligible doublon fraction, attractive interactions( a s <
0) favor initial states with a discernible doublonfraction, as illustrated in Fig. S1(a). Creating initialstates with a large doublon fraction, however, also in-creases the number of holes, since a pair of neighboringsinglons is converted into a doublon and a hole. In or-der to generate initial states with an appreciable frac-tion of atoms on doubly-occupied sites [ n d = 0 . a s = − a . Thisvalue was chosen in order to maximize the doublon frac-tion, while at the same time minimizing the cloud ra-dius and therefore maximizing the density. As shownin Fig. S1(b), the cloud radius is the same as the one,where the scattering length was set to 25 a during load-ing; with U/ (12 J ) = 0 .
12 in the (8 E rx , E r ⊥ , E r ⊥ ) deeplattice. The measured cloud radius [Fig. S1(b)] is inagreement with the values obtained in Ref. [63], where,for our parameters, the central density of the cloud wasshown to be (cid:104) ˆ n i (cid:105) (cid:46) .
9. For initial states without dou-blons ( n d < . a s = 140 a . S3. Creating a flat potential for the expansion
At the end of the initial state preparation, the atomiccloud is confined by three dipole traps: one along each ofthe horizontal directions, x and y ( x is the longitudinaldirection of the tubes), and one along the vertical z direc-tion. The vertical dipole trap has a Gaussian beam waistof 150 µ m. The horizontal dipole traps are elliptical withwaists of 30 µ m in the vertical and waists of 300 µ m in thehorizontal direction. The optical lattices along all spatialaxes have beam waists of 150 µ m (as the vertical dipoletrap) and are blue detuned, providing an anti-confiningpotential. A flat potential along the x direction duringthe expansion of the cloud can be generated by choos-ing the strength of the vertical dipole trap such that itcompensates the anti-confinement of the optical lattices.As the horizontal dipole traps have different beam ge-ometries, they cannot be used to compensate the anti-confinement, and therefore, they are switched off duringthe expansion measurements. Creating a flat potentialrequires optimizing both the z dipole beam alignmentand its strength to maximize the in-situ cloud size af-ter a long evolution time (see also Ref. [19]). However,a completely flat potential cannot be implemented usingthis method, as the harmonic contributions of the overallpotential can only be canceled within a certain area. S4. Singlon- and doublon-resolved analysis
Doublons can be removed from the lattice using an ad-ditional laser pulse, which is blue-detuned by an amount∆ from the imaging transition of K ( | F = 9 / , m F = − / (cid:105) → | F (cid:48) = 11 / , m F (cid:48) = − / (cid:105) ). As a result, atoms Pulse duration (ms) N s + N d FIG. S2.
Atom number decay in the presence of anadditional near-resonant laser pulse:
Total atom number N s + N d as a function of the pulse duration. The data wastaken for U = 20 J after an expansion time in the lattice of t = 40 τ at a detuning of the pulse of ∆ = 364 MHz. Solidlines show the fit using the sum of two exponentials. on doubly-occupied sites are lost due to light-assisted col-lisions. As indicated by the semi-log-plot in Fig. S2, thetotal atom number shows a bimodal decay with well-separated time scales, which we characterize by fittinga sum of two exponential functions to the atom num-ber N s ( t ) + N d ( t ) = N s ( t = 0) exp (cid:0) − t/τ s (cid:1) + N d ( t =0) exp (cid:0) − t/τ d (cid:1) . We extract a fast decay with a lifetime of τ d = 40(10) µ s, which is attributed to the loss of atoms ondoubly-occupied sites due to light-assisted collisions andan additional slow decay with a lifetime of τ s = 12(1) ms,which results from the loss of atoms on singly-occupiedsites due to off-resonant photon scattering. For all exper-iments, the duration of the light pulse was set to 150 µ s.The detuning ∆ of the pulse with respect to resonantimaging light depends on the magnetic field used to setthe final interaction strength with the Feshbach reso-nance and varies between ∆( U = 5 J ) = 296 MHz and∆( U = 20 J ) = 364 MHz. We found that the differentdetunings have only a negligible effect on the doublonlifetimes. By subtracting the optical density of succes-sive pictures with and without light pulse, doublon- andsinglon-resolved dynamics can be analyzed individually. S5. Estimation of tubes with few doublons
Absorption imaging along the z direction results in av-eraging over many individual realizations of 1D systemswith different density distributions. Therefore, the in-crease of N cd /N cs in Fig. 2(d) of the main text could inprinciple be caused by 1D systems that contain singlonsbut only a very small amount of doublons. In these sys-tems, the expansion would decrease the singlon numberin the center N cs , without mediating doublon dynamics.3 R d ( d ) U =5 J U =7.5 JU =10 J U =15 J (a)(b) R d ( d ) Expansion time t ( ¿ )Expansion time t ( ¿ ) U =20 J FIG. S3.
Time traces of R d with linear fits : (a) Timetraces of the HWHM of the doublon cloud R d for differentinteractions. Solid lines indicate a linear fit to the data. (b)Time trace of R d for U = 20 J with a linear fit resulting in χ = 0 .
61 (left) and a fit of a constant function (dashedline) resulting in χ = 0 .
84 (right).
In order to estimate the contribution of 1D systems withnegligible doublon fraction to the ratio N cd /N cs , we canapproximate the initial shape of the cloud in the latticewith an ellipsoid, where the ratio of the principal axesis set by the trap frequencies. Using an Abel transform,the full 3D density distribution of singlons and doublonscan be reconstructed. The width of the singlon and dou-blon distribution along the x direction thereby translatesinto a width along the z direction by implying the sym-metry of the ellipsoid. As a result, the extent of thecloud along the z direction amounts to about 30 indi-vidual planes, where only the four outermost planes areprimarily occupied with singlons. Assuming that the dy-namics and therefore the change in the ratio N cd /N cs issolely governed by singlons expanding in these outermost1D systems we obtain a conservative upper bound of N cd /N cs (cid:46) S6. Fits to the time traces of R d In order to quantify the dynamics of the width of thedoublon cloud in Fig. 2, we fit the time traces of R d with a linear function f ( t ) = at + b . The fits to the data areshown in Fig. S3(a) for increasing interactions and addi-tionally in the left panel of Fig. S3(b) for the strongestinteraction of U = 20 J . While the fits indicate an in-crease of R d for weak interactions of U = 5 J , strongerinteractions result in a slower spreading of the doubloncloud and for U = 20 J , the fit clearly indicates a con-traction of the doublon cloud as shown in the left panelof Fig. S3(b). In order to test the goodness of the linearfit to the doublon cloud size for strong interactions of U = 20 J , we compare it to fitting a constant function tothe doublon cloud size in the right panel of Fig. S3(b). Acommon quantitative estimate for the goodness of a fit is χ , which is defined as χ = ν − N (cid:88) i =1 (cid:18) R d ( t i ) − g ( t i ) σ i (cid:19) , (S1)where R d ( t i ) is the doublon cloud size at time t i , calcu-lated from the mean of four data points, σ i is the corre-sponding standard deviation for each R d ( t i ), g ( t i ) is thevalue of the fit function at time t i and ν is the differ-ence between the total number N of discrete points t i ( N = 17 for all time traces of the doublon cloud size)and the number of fitting parameters. For the linear fit,we obtain χ = 0 .
61, whereas the fit with the constantfunction yields a slightly larger χ = 0 .
84. This χ analysis indicates that the decreasing linear function de-scribes the data better than a constant function, evenafter accounting for the increased number of fit parame-ters. Using the resulting fit parameters, we calculate therelative change in doublon cloud size ∆ R d = a/b · τ for all interactions and show the results in the inset ofFig. 3 in the main text. The error of the linear fittingparameters indicates one sigma confidence intervals andthe errorbars of the relative doublon cloud size are cal-culated by using Gaussian error propagation. S7. Extracting radial velocities
The size of the cloud can be characterized by mea-suring the second moment r , which is defined as r = (cid:80) l ρ l ( l − l c ) · . d , where ρ l is the normalized opticaldensity at pixel l with distance | l − l c | from the centralpixel l c of the 1D integrated line densities and the factor8 . d . The second moment is less affected by detailsof the density distribution than R d and R s , since it takesinto account the density of the whole cloud and not onlythe height of the cloud at a chosen width. This, however,leaves it more susceptible to noise, compromising its ap-plicability in the analysis of the doublon clouds (Fig. 2 ofthe main text). The second moment r is extracted bysubtracting the background from the raw integrated linedensities and then summing over the integrated line den-sities of the cloud from the center outwards, until r sat-urates. The cloud size at every time step is averaged over4 r ( d ) Expansion time t ( ¿ ) r ( d ) (a)(b) U =5 JU =20 J FIG. S4.
Second moments r : Exemplary time traces ofthe cloud size r for (a) U = 5 J and for (b) U = 20 J . Solidlines are fits to the time traces to extract the radial velocities. five measurements and error bars denote the standard er-ror of the mean. Exemplary time traces of r = √ r areshown in Fig. S4(a) for U = 5 J and in Fig. S4(b) for U = 20 J . From the time dependence of r , we extractthe radial velocity v r by fitting √ r = (cid:112) r + v r t , where r is the second moment of the cloud at t = 0. We do notuse ˜ r as defined in Eq. (S4) in the following Sec. S8.3,since this would require a sufficiently accurate knowledgeof r . In the non-interacting case, the radial velocity of aninitially localized particle yields v r = √ d/τ (in agree-ment with our experimental results for U = 0) and arean average in quasi-momentum space over the group ve-locities weighted with the momentum-distribution func-tion [32, 33]. S8. tDMRG simulations
1. Method
The non-equilibrium expansion is studied by meansof the time-dependent density-matrix renormalizationgroup (tDMRG) method (see Refs. [64–67] for details) us-ing a Trotter-Suzuki decomposition scheme for the timepropagation. The results presented in Fig. 3 of the maintext were obtained with a second-order Trotter-Suzukidecomposition, up to M = 6000 DMRG states, a dis-carded weight of 10 − and δt = 0 . τ . For the datashown in Figs. S7 and S8, δt = 0 . τ was used. We en-sure the numerical accuracy of the results by varying thetime step, the discarded weight during the time evolu-tion and the maximum number of states. Some of ourtDMRG implementations make use of the tensor librarydeveloped in [71].The initial states are prepared on a lattice of L sites as a product state, i.e., | ψ (cid:105) = (cid:89) i = i L ,...,i R (cid:16) ˆ c † i ↑ (cid:17) n i ↑ (cid:16) ˆ c † i ↓ (cid:17) n i ↓ | (cid:105) , (S2)with i L = L/ − L init / i R = L/ L init /
2, where L init is the size of the region with a nonzero density at t = 0and n iσ ∈ { , } are specially chosen integers. Thus, eachsite inside the initially confined region is occupied by ex-actly one singlon, doublon, or holon (empty site). Theexperimental results are averages over many 1D tubesin each measurement of the cloud and, as a consequence,various singlon/doublon/holon configurations are probedsimultaneously. In order to account for this, the resultspresented in Fig. 3 of the main text are averaged over 120simulations, each starting from different (random) distri-bution of particle configurations. The latter ensures anapproximately uniform density in the initially ( t = 0)confined region. Note that the results contain an errordue to the sampling over a subset of all possible configu-rations. This small error is comparable in magnitude toour numerical accuracy at the maximal considered time(order of a few percent). In addition, each sample is pre-pared with a random distribution of spin-up and spin-down fermions under the constraint of N ↓ = N ↑ . Due tothe entanglement increase in tDMRG simulations [67], itis unfeasible to simulate the full distribution of 1D tubesrealized in the experiment for realistic particle numbersand time scales.
2. Initial states with doublons
We keep the total number of doublons fixed to N d / (cid:80) i (cid:104) ˆ n i ↑ ˆ n i ↓ (cid:105) = 4 and set the singlon-to-doublon ratio to2 N s /N d = 3, where N s = N − N d with N = (cid:80) iσ (cid:104) ˆ n iσ (cid:105) .Note that N d denotes the number of atoms on doubly-occupied sites, hence, the number of doublons is N d / L = 60 sites and we only consider timesbefore reflections off the boundaries occur.As is clearly visible in Fig. 3 of the main text, the den-sity of the initial state highly influences the quantum-distillation dynamics. Here, we show how additional de-tails of the initial cloud affect the achievable doubloncontraction and the transient expansion dynamics. Inorder to account for typical aspects of the experiment,we consider four types of initial states (see also Fig. S5): A : A box trap with 12 singlons and 4 doublons and aninitial box size of L init = 16, B : A box trap with 12 singlons, 4 doublons, and 4holons and an initial box size of L init = 20, C : A box trap with 12 singlons, 4 doublons, and 8holons and an initial box size of L init = 24, D : case A surrounded by singlon wings, i.e., two addi-tional singlons placed on the left and the right sideof the configurations used in A ,5 ABCD ̂ n i ®® ̂ n i ®® ̂ n i ®® ̂ n i ®® iiii FIG. S5.
Schematics of the different initial states withdoublons : A : Box trap with singlons and doublons, B and C : Box trap containing doublons, singlons and holons, D :singlons and doublons as in case A but surrounded by sin-glon wings. The shaded area depicts the region in which werandomly distribute doublons, singlons, or holons, keeping N d / N s = 12 fixed. Outside those regions, we onlyallow for singlons (case D ) or no atoms. Case A corresponds to the box trap with solely sin-glons and doublons. This idealized setup was previouslyinvestigated in theoretical studies [42, 44]. Cases B and C also include 4 or 8 holons, respectively, and thus L init increases as we go from A to C . At the same time, theaverage density goes down. We choose to keep N d / R d in Fig. 3 of the main text. The presence ofholons can be viewed as imperfections resulting duringthe loading process and initial-state preparation. Case D accounts for the inhomogeneous shape of the initialfermionic cloud resulting from the trapping potential.Complementary to the main text, where we show theHWHM R d , here, we discuss the doublon cloud radius de-fined as r d ( t ) = (cid:113) (1 /N d ( t )) (cid:80) i n di ( t )( i − i ) , with i thecenter of mass (here, i = L/ .
5) and n di = (cid:104) ˆ n i ↑ ˆ n i ↓ (cid:105) is averaged over all random configurations. In the mainpanel of Fig. S6, we present the relative change of thecloud radius ∆ r d ( t ) = r d ( t ) /r d (0) − U/J = 20. We stress the main observations:(i) The time traces of ∆ r d behave qualitatively very sim-ilar to the time traces of ∆ R d (relative change in theHWHM shown in the inset of Fig. S6). (ii) The pres-ence of holons, cases B and C , reduces the achievableminimal cloud size on the time scales of the quantumdistillation [44]. (iii) The presence of additional singlonsat the edges of the cloud (case D ) allows doublons tomove outwards on the singlon wings in the early stage ofthe expansion. As a consequence, the presence of singlonwings delays the contraction of doublons. Expansion time t ( ¿ ) −0.5−0.4−0.3−0.2−0.10.00.1 ¢ r d A : n =1.25 B : n =1.00 C : n =0.83 D : A with singlon wings 0 5 10 t ( ¿ ) ¢ R d FIG. S6.
Doublon cloud radius : Time evolution of ∆ r d for U = 20 J and various initial configurations. Inset: Timetraces of ∆ R d (relative change in the HWHM).
3. Initial states without doublons
The initial states without doublons are prepared byenforcing n i ↑ + n i ↓ = 1 in Eq. (S2) on a lattice of L =100 sites with N = 10 particles placed in the center ofthe system. Note that in this case, we keep the particlenumber fixed and do not consider holons in the initialstate, while our results are averaged over all possible spinconfigurations with N ↑ = N ↓ .In order to extract the velocities, we first calculate ˜ r ( t )via: r ( t ) = 1 N L (cid:88) i =1 (cid:104) ˆ n i (cid:105) ( i − i ) , (S3)˜ r ( t ) = (cid:112) r ( t ) − r (0) , (S4)where (cid:104) ˆ n i (cid:105) is the density at site i averaged over all initialspin configurations. The procedure to extract the asymp-totic velocities is illustrated in Fig. S7 where the shadedregion indicates the fitting window. Clearly, ˜ r is essen-tially linear in time and in the fitting window the velocityhas settled to a constant value, as shown in Fig. S7(b).The different initial states can be discriminated by thenumber of domain walls that is the number of times anup-spin particle is next to a down-spin particle. For N = 10 particles in the initial state, the minimum num-ber of domain walls is one and the maximum number isnine. We calculate the velocities for different numbersof domain walls in the initial state as before by first av-eraging over all configurations with the same number ofdomain walls and then extract v r . In Fig. S8, we plotthe velocities for different numbers of domain walls ver-sus the interaction energy at time tJ = 8 as diamonds(one to nine domain walls in the initial state from left to6 ̃ r ( d ) Expansion time t ( ¿ ) d ̃ r / d t ( d / ¿ ) (a)(b) FIG. S7.
Extraction of the radial velocity v r from thetDMRG data : (a) ˜ r as defined in Eq. (S4) as a function oftime indicated as a red solid line and linear fit to the curveindicated by a dashed black line. (b) Time derivative of ˜ r to determine the time when the transient behavior is finishedand ˜ r is approximately linear in time. The shaded regionrepresents the time interval where the linear fit is done. right) [33]. The interaction energy is defined as: E int = U L (cid:88) i =1 (cid:104) ˆ n i ↑ ˆ n i ↓ (cid:105) . (S5)We also carried out a simulation of the Bose-Hubbardmodel where we start with a region of 10 singly-occupiedsites in the middle of an otherwise empty lattice [19, 33].For the Bose-Hubbard model, this is the unique productstate with one boson per site and doublons and highersite occupancies can be generated dynamically in everysite. The interaction energy for the Bose-Hubbard modelis defined as: E int = U L (cid:88) i =1 (cid:104) ˆ n i (ˆ n i − (cid:105) . (S6)The data for the Bose-Hubbard model is shown as cir-cles in Fig. S8. The interaction quench from U/J = ∞ to U/J < ∞ causes the dynamical formation of dou-blons (there were none in the initial state). The trapopening induces a decrease of E int towards the asymp-totic value. We do not reach the asymptotic regime in our simulations, but we choose a time large enough tocapture most of the decay of E int . The results in Fig. S8suggest that for large U/J , the asymptotic radial velocityis indeed primarily a function of the interaction energythat is generated due to the interaction quantum quenchover the first tunneling times [69, 70]. There is a notice-able and expected additional U -dependence (see the insetwhere we plot v r versus E int /U ). This results from (i)the fact that doublons are only well-conserved objects for U (cid:29) W and (ii) that doublons expand with a nonzero ve-locity for finite U ∼ W . The argument of Ref. [57], whichexplains the expansion velocities at large U/J , assumesimmobile doublons (and higher site occupancies) on therelevant time scales, which is only correct for U (cid:29) W .The interaction energy is only a proxy for the actualheavy objects involved in the dynamics, in particular,since E int still undergoes a slow decrease beyond thetimes reached in the simulations. A more rigorous ar-gument is to relate v r to the overlap of the initial statewith bound states (see [68] for the two-body case) in theintegrable 1D FHM model, in extension of the approachtaken in [35]. E int ( tJ =8)/ J v r ( d / ¿ ) U =20 JU =10 JU =5 J E int ( tJ =8)/ U v r ( d / ¿ ) FIG. S8.
Radial velocity as a function of the in-teraction energy for different interaction strengths :The three colors correspond to interaction strengths
U/J =5 , ,
20. Diamonds correspond to fermions, circles corre-spond to bosons. Diamonds of the same color correspondto different numbers of domain walls in the initial state (onlyone domain wall to nine domain walls from left to right). Solidlines are quadratic fits to the data. The inset shows the samedata but versus E int /U/U