Twofold correlation spreading in a strongly correlated lattice Bose gas

We study the spreading of correlations in the Bose-Hubbard chain, using the time-dependent matrix-product state approach. In both the superfluid and the Mott-insulator phases, we find that the time-dependent correlation functions generally display a universal twofold cone structure characterized by two distinct velocities. The latter are related to different microscopic properties of the system and provide useful information on the excitation spectrum. The twofold spreading of correlations has profound implications on experimental observations that are discussed.

In this supplemental material, we give more details about several points discussed in the main paper. In Sec. S1, we discuss the time-dependent matrix-product state (t-MPS) simulations and the choice of the numerical parameters to ensure the convergence of the results in all the regimes of the Bose-Hubbard model considered in the main paper. In Sec. S2, we present t-MPS results for the spreading of the one-body correlation function G 1 in the mean-field superfluid (SF) regime. Section S3 briefly outlines the mapping from the 1D Bose-Hubbard model to the Lieb-Liniger model and gives the correspondance of the parameters. Finally, in Sec. S4 we discuss the strong-coupling expansion of the correlation function G 2 for unit filling, n = 1, and discuss the suppression of its twofold structure deep in the Mott insulator (MI) phase.

S1. TIME-DEPENDENT MATRIX-PRODUCT STATE SIMULATIONS
The numerical results reported in the main paper are all obtained using the time-dependent density-matrix renormalization group approach (DMRG) with the matrix-product state representation (t-MPS approach) [1][2][3]. It yields numerically-exact results on both equilibrium and out-of-equilibrium properties of low dimensional lattice models. The approach resorts on the Schmidt expansion of the many-body wave function and permits to reduce the Hilbert space to a finite, relevant subset, provided the entanglement entropy remains sufficiently small. Owing to the area law [4,5], it is optimal for 1D lattice models with a finite local Hilbert space in gapped phases, the entanglement of which remains finite in the thermodynamic limit. It also applies to gapless phases, although with more stringent numerical parameters (high-filling cut-off and the bond dimension). To validate the accuracy of our results in all phases of the BH model, a systematic study of the effect of these parameters has been performed.
Truncation of the local Hilbert space.-For the BH model considered in this work, the local Hilbert space is spanned by the Fock basis of number states, |n R , where n R ∈ N, which is infinite. However, the probability distribution of the latticesite occupation n R decays faster than exponentially in both the SF and MI phases. Accurate results can thus be obtained by cutting off the local Hilbert space to some value n max . It is important to note that, in some cases, the value of n max needs to be significantly much larger than the average filling n and its fluctuations. This observation is consistent with analyses of truncated Bose-Hubbard models using quantum Monte Carlo simulations [6].
The SF mean-field regime, which corresponds to a high filling factorn and the gapless dispersion relation, has the most binding criteria. We found that a good estimator for n max is given by the condition 1 − nmax n=0 P (n) 10 −2 , where P (n) is the probability that n bosons occupy a given lattice site. In the SF mean-field regime, the probability distribution is nearly Poissonian, P (n) n n e −n /n!. For instance, for the filling factor n = 5 used for the data of Fig. 2, it yields n max 12. For the strongly correlated SF regime at n = 1 considered for Fig. 3(a), the density fluctuations are significantly suppressed and using the same condition as previously leads to n max = 5. For the MI phase atn = 1 and moderate values of U/J (15 ≥ U/J ≥ u c ) considered for Fig. 3(b), we kept n max = 5. Deep in the MI phase (U/J ≥ 15), truncating the local Hilbert space to n max = 2, as used for Fig. 3(c) turns out to be sufficient. Finally, the strongly interacting SF regime is the easiest case from a numerical point of view. Owing to the low filling factorn < 1 and the large value of the interaction parameter U/J, the above condition also yields n max = 2, as used for Fig. 4. In all cases, we have checked that the numerics are converged for these values of n max .
Bond dimension.-Within the MPS approach, the many-body state for a M -site lattice is represented in the tensor network form where n j spans a local Hilbert space basis. For the BH model, it corresponds to a Fock basis truncated at n max . For each value of n j , the quantity A nj [j] is a χ j−1 × χ j matrix, where χ j is the rank associated to the Schmidt matrix when applying the j-th singular value decomposition [2]. The bond dimension χ is defined as the maximum rank, Note that for open-boundary conditions, the quantities A n1 [1] and A n M [M ] are actually a row vector and a column vector, respectively, i.e. χ 0 = χ M = 1.
In the numerics, the maximum value of χ is chosen sufficiently large so that the truncation does not affect the results. In practice, the calculations are run for several values of χ up to convergence of the correlation function G 1 (R, t) or G 2 (R, t). The required value of χ significantly depends on the regime and on the observable. In the following, we give the values used for the final results presented in the paper.
For the SF mean-field regime [Figs. 2(a) and S1], we used the values χ = 300 and χ = 450 for the G 2 and G 1 functions, respectively. The bond dimension used for G 1 is higher than the one for G 2 due to the long-range phase correlations already present at equilibrium. For the SF strongly correlated regime at n = 1 [ Fig. 3(a)], we used χ = 300 for both correlation functions. A similar value of χ was considered for moderate values of U/J in the MI phase atn = 1 [ Fig. 3(b)]. Deep in the MI phase [ Fig. 3(b)], the bond dimension can be significantly decreased and we consider χ = 100. Finally, in the SF strongly interacting regime at U/J = 50, we found that the value χ = 100 is enough.

S2. ONE-BODY CORRELATION FUNCTION G1(R, t) IN THE MEAN-FIELD REGIME
In the analysis of the SF mean-field regime reported in the main paper, we focused on the two-body correlation function G 2 (R, t). We have also studied the one-body correlation G 1 (R, t) using the same t-MPS simulations. We found that the dynamics of the G 1 function shows a spike-like structure, similar to that found for the G 2 function. The values of the correlation edge (V CE ) and maxima (V m ) velocities agree with those found for the G 2 function within less than 10%. Figure S1 shows an example, for the quench from (U n/J) 0 = 1 to U n/J = 0.5, andn = 5. The fits to the correlation edge and to the maxima yield the velocities V CE = (4.4 ± 0.3) J/ and V m = (3.3 ± 0.2) J/ , in excellent agreement with the corresponding values found from the dynamics of the G 2 function, see Fig. 2(b).
The agreement between the spreading velocities for different correlation functions was found in all regimes, see for instance Figs. 3(d1) and (d2). It is consistent with the prediction that these velocities are characteristic of the excitation spectrum and not on the details of the correlation function [7]. Note, however, that the full space-time dependence of the signal depends on the correlation function. In general, we found that the signal for G 1 is less sharp than for G 2 . This may be attributed to the long-range phase correlations present in the initial state, which blur the correlation function [8]. Figure S1: Spreading of the one-body correlation function G1(R, t) for a global quench in the SF mean-field regime from (U/J)0 = 0.2 to U/J = 0.1 andn = 5. The solid-green and dashed-blue lines are fits to the CE and maxima, respectively.

S3. MAPPING ON THE 1D LIEB-LINIGER MODEL
In the long-wave length regime, the lattice discretization of the Bose-Hubbard (BH) may be disregarded. The BH model then maps onto the continuous-space Lieb-Liniger (LL) model, It describes a one-dimensional gas of N bosons of mass m with contact interactions, characterized by the interaction strength c > 0. The correspondance between the parameters of the BH and LL models is found by discretizing the LL model, Eq. (S2), on the length scale defined by the lattice spacing a. It yields J = 2 /2ma 2 and U = 2 c/ma. The density of the LL model is ρ ≡ N/L = n/a, where n is the number of bosons per lattice site (filling) and L is the system size.
The LL Hamiltonian is exactly solvable by Bethe ansatz [9,10]. All the thermodynamic quantities at zero temperature can be written as universal functions of the Lieb-Liniger parameter γ = c/ρ and the dimensionless quantity e(γ) = E 0 /N n 2 , where E 0 is the ground state energy. For instance, the macroscopic sound velocity [10] reads as v s ≡ L mρ Using the small γ expansion, e(γ) = γ 1 − (4/3π) √ γ , one then finds valid in the weakly-interacting regime, γ 1. Finally, using the correspondance between the parameters of the BH and LL models, one finds and γ = U/2Jn.

S4. TWO-BODY CORRELATION FUNCTION G2(R, t) IN THE MOTT-INSULATING PHASE
In order to explain the suppression of the twofold structure for the two-body correlations deep in the Mott insulator phase (MI; U J and n = 1), we compute the function G 2 (R, t), working along the lines of Ref. [11]. Considering the manifold of doublon-holon pairs and mapping the resulting Hamiltonian into a fermionic one, the two-body correlation function may be written as with and the excitation spectrum is 2E k [U − 2J(2n+1) cos(k)] 2 + 16J 2n (n+1) sin 2 (k), see Eq. (3).
Quench deep into the Mott insulator phase.-For a quench, very deep in the MI phase, U J, the second right-hand-side term in Eq. (S6) is much smaller than the first one and the former can be neglected. Using Eq. (S7), it yields explicitly for G 2 (R, t) −2|g 2 (R, t)| 2 , Moreover, the excitation spectrum may be expanded in powers of J/U . Up to first-order, it yields 2E k U −2J(2n+1) cos(k).
The gap term e iU t can then be factorized in the two terms under the integral in Eq. (S9) and disappears due to the square modulus. Introducing the effective excitation spectrum 2Ẽ k = −2J(2n + 1) cos(k), we then find G 2 −2|g 2 (R, t)| 2 with The integral may be evaluated using the stationary phase approximation. In the infinite time and distance limit along the line R/t = cst, the integral in Eq. (S10) is dominated by the momentum contributions with a stationary phase (sp), i.e. ∂ k (2Ẽ k t ± kR) = 0 or, equivalently, 2Ṽ g (k sp ) = ±R/t whereṼ g = ∂ kẼk is the group velocity of the effective excitation spectrum. Since the latter is upper bounded by the valueṼ * g = max(Ṽ g ) = J(2n + 1), it has a solution only for R/t < 2Ṽ * g . We then find g 2 (R, t) ∼ J UṼ g (k sp ) |∂ 2 kẼ ksp |t 1/2 cos 2Ẽ ksp t − k sp R + σ π 4 + i sin 2Ẽ ksp t − k sp R + σ π 4 . (S11)