Twofold correlation spreading in a strongly correlated lattice Bose gas

TL;DR

This paper investigates how correlations spread in a lattice Bose gas, revealing a universal twofold cone structure with two distinct velocities in both superfluid and Mott-insulator phases, providing insights into the excitation spectrum.

Contribution

It demonstrates the universal twofold correlation spreading in a strongly correlated Bose-Hubbard chain using advanced numerical methods.

Findings

Identification of a twofold cone structure in correlation spreading

Distinct velocities related to microscopic properties

Implications for experimental observations

Abstract

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.

Equations

$$\hat{H}=-J\sum _R\left({\hat{a}}_R^{†}{\hat{a}}_{R+1}+h.c.\right)+\frac{U}{2}\sum _R{\hat{n}}_R({\hat{n}}_R-1),$$

$$E_k\simeq \sqrt{{\epsilon }_k\left({\epsilon }_k+2\bar{n}U\right)},$$

$$2E_k\simeq \sqrt{{\left[U-2J(2\bar{n}+1)\cos (k)\right]}^2+16J^2\bar{n}(\bar{n}+1){\sin }^2(k)}.$$

$$g_{2}(R,t)\propto\frac{J}{U}\frac{R}{t}\int_{-\pi}^{+\pi}\frac{\textrm{d}k}{2\pi}\Big{\{}\textrm{e}^{i\left(2E_{k}t+kR\right)}+\textrm{e}^{i\left(2E_{k}t-kR\right)}\Big{\}}.$$

$$|\Psi ⟩=\sum _{n_1,n_2,\dots n_M}{A}^{n_1}[1]A^{n_2}[2]\dots A^{n_M}[M]|n_1,n_2,\dots ,n_M⟩,$$

$$\hat{H}=\frac{{\hslash }^2}{2m}\left[-\sum _{i=1}^N\frac{{\partial }^2}{\partial x_i^2}+c\sum _{i\ne j}\delta (x_i-x_j)\right].$$

$$v_{\text{s}}\equiv \sqrt{\frac{L}{m\rho }{\frac{{\partial }^2{E}_0}{\partial L^2}|}_{N,S}}=\frac{\hslash \rho }{m}\sqrt{3e(\gamma )-2\gamma e^{\prime }(\gamma )+\frac{1}{2}{\gamma }^2{e}^{\prime \prime }(\gamma )}.$$

$$v_{\textrm{s}}=\dfrac{\hbar\rho}{m}\sqrt{\gamma}\big{(}1-\sqrt{\gamma}/4\pi\big{)},$$

$$V_{\textrm{s}}\equiv v_{\textrm{s}}/a=\frac{2J{\overline{n}}}{\hbar}\sqrt{\gamma}\big{(}1-\sqrt{\gamma}/4\pi\big{)}$$

$$G_{2}(R,t)\simeq-2\big{(}|g_{2}(R,t)|^{2}+|\bar{g}_{2}(R,t)|^{2}\big{)},$$

$$\displaystyle g_{2}(R,t)\sim\frac{J}{U}\frac{R}{t}\int_{-\pi}^{+\pi}\frac{\textrm{d}k}{2\pi}\Big{\{}\textrm{e}^{i\left(2E_{k}t+kR\right)}+\textrm{e}^{i\left(2E_{k}t-kR\right)}\Big{\}},$$

$$\displaystyle\bar{g}_{2}(R,t)\sim\left(\frac{J}{U}\right)^{2}\int_{-\pi}^{+\pi}\frac{\textrm{d}k}{2\pi}\sin^{2}(k)\Big{\{}\textrm{e}^{i\left(2E_{k}t-kR\right)}+\textrm{e}^{-i\left(2E_{k}t+kR\right)}\Big{\}}$$

$$G_{2}(R,t)\sim-2\left(\frac{J}{U}\right)^{2}\left(\frac{R}{t}\right)^{2}\left|\int_{-\pi}^{\pi}\frac{\textrm{d}k}{2\pi}\Big{\{}\textrm{e}^{i\left(2E_{k}t+kR\right)}+\textrm{e}^{i\left(2E_{k}t-kR\right)}\Big{\}}\right|^{2}$$

$$g_{2}(R,t)\sim\frac{J}{U}\frac{R}{t}\int_{-\pi}^{\pi}\frac{\textrm{d}k}{2\pi}\Big{\{}\textrm{e}^{i\left(2\tilde{E}_{k}t+kR\right)}+\textrm{e}^{i\left(2\tilde{E}_{k}t-kR\right)}\Big{\}}.$$

$$g_2(R,t)\sim \frac{J}{U}\frac{{\tilde{V}}_{\mathrm{g}}(k_{\mathrm{sp}})}{{\left(|{\partial }_k^2{\tilde{E}}_{k_{\mathrm{sp}}}|t\right)}^{1/2}}\left[\cos \left(2{\tilde{E}}_{k_{\mathrm{sp}}}t-k_{\mathrm{sp}}R+\sigma \frac{\pi }{4}\right)+i\sin \left(2{\tilde{E}}_{k_{\mathrm{sp}}}t-k_{\mathrm{sp}}R+\sigma \frac{\pi }{4}\right)\right].$$

$${\bar{g}}_2(R,t)\sim {\left(\frac{J}{U}\right)}^2\frac{{\sin }^2(k_{\mathrm{sp}})}{{\left(|{\partial }_k^2{E}_{k_{\mathrm{sp}}}|t\right)}^{1/2}}\cos \left(2E_{k_{\mathrm{sp}}}t-k_{\mathrm{sp}}R+{\sigma }^{\prime }\frac{\pi }{4}\right)$$

Full text

In the last decades, simultaneous progress of the many-body quantum theory and the experimental control of quantum matter in condensed matter and atomic, molecular, and optical physics has given dramatic momentum to the understanding of the out-of-equilibrium dynamics of correlated quantum systems polkovnikov2011 ; eisert2015 ; *pekola2015; langen2015 ; lewenstein2007 ; *bloch2008; NaturePhysicsInsight2012bloch ; *NaturePhysicsInsight2012blatt; *NaturePhysicsInsight2012aspuru-guzik; *NaturePhysicsInsight2012houck. The spreading of quantum correlations governs many fundamental phenomena, including the propagation of information and entanglement, thermalization, and the area laws for entanglement entropy. For lattice systems with local interactions, the existence of Lieb-Robinson (LR) bounds implies the emergence of a causal light cone beyond which the correlations are exponentially suppressed  lieb1972 ; bravyi2006 ; hastings2006 . So far, light-cone-like spreading of correlations has been reported in short-range interacting models barmettler2012 ; cheneau2012 ; carleo2014 ; manmana2009 as well as long-range models jurcevic2014 ; richerme2014 ; hauke2013 ; eisert2013 ; cevolani2015 ; schachenmayer2015b ; buyskikh2016 ; cevolani2016 ; frerot2018 ; cevolani2018 where weaker LR bounds exist hastings2006 ; foss-feig2015 . However, many questions remain open. For instance, it is still debated whether a non-linear cone emerges in generic long-range systems, for which different results point towards either super-ballistic, ballistic or sub-ballistic spreading. It was recently proposed that these apparently conflicting results can be reconciled by the coexistence of several signals governed by different scaling laws cevolani2018 . This behavior may be related to the non-linearity of the quasiparticle excitation spectrum, and may also appear in systems with short-range interactions. In the later case, it is expected that both the signals spread ballistically but with different velocities. However, this picture relies on meanfield theory, which ignores potentially important dynamical effects, such as quasiparticle collisions and finite lifetime.

In this work, using an exact many-body approach beyond meanfield theory, we demonstrate the emergence of a universal twofold dynamics in the spreading of correlations for a generic short-range, strongly correlated quantum model. Specifically, we consider the one-dimensional Bose-Hubbard model and use time-dependent tensor network techniques based on matrix product states. Spanning the phase diagram, we almost always find a twofold structure of the space-time correlation pattern, characterized by two distinct velocities, essentially irrespective of the correlation function. Exceptions, discussed below, only appear for particular cases. In the superfluid meanfield regime and in the Mott insulator phase, the two spreading velocities are readily interpreted from the properties of the corresponding excitation spectra, which are known. In the strongly correlated superfluid regime, only the sound velocity is known. There, our results show beyond Luttinger liquid behavior and provide useful information about the excitation spectrum beyond the phonon branch. The emergence of a universal twofold spreading of correlations has profound implications on experimental observations, which we discuss, including with a view towards extensions to long-range systems.

Model and approach.—

The Hamiltonian of the one-dimensional (1D) Bose-Hubbard (BH) model, considered throughout this work, reads as

[TABLE]

where $\hat{a}_{R}$ and $\hat{a}_{R}^{\dagger}$ are the bosonic annihilation and creation operators on site $R$, $\hat{n}_{R}=\hat{a}^{\dagger}_{R}\hat{a}_{R}$ is the occupation number (filling), $J$ is the hopping amplitude, $U>0$ is the repulsive on-site interaction energy, and the lattice spacing is fixed to unity ($R\in\mathbb{Z}$). At equilibrium and zero-temperature, the phase diagram of the 1D BH model is well known  sachdev2001 ; cazalilla2011 , and sketched on Fig. 1(a). It comprises a superfluid (SF) and a Mott insulator (MI) phase, determined by the competition of the hopping, the interactions, and the average filling ${\overline{n}}$ (or, equivalently, the chemical potential $\mu$). For commensurate filling, ${\overline{n}}\in\mathbb{N}^{*}$ the SF-MI (Mott-$U$) transition is of the Berezinskii-Kosterlitz-Thouless type, at the critical value $u_{\textrm{c}}\simeq 3.3$ for unit filling (${\overline{n}}=1$) in 1D kuhner2000 ; kashurnikov1996exact ; ejima2011 ; rombouts2006 . For incommensurate filling, the Bose gas is a SF for any value of $U/J$. The commensurate-incommensurate (Mott-$\delta$) transition, of the meanfield type, is then driven by doping when ${\overline{n}}$ approaches a positive integer value for sufficiently strong interactions.

We study the out-of-equilibrium dynamics of the BH model by applying a sudden global quench calabrese2006 ; barmettler2012 ; kollath2007 ; moeckel2008 ; manmana2009 ; roux2010 ; navez2010 ; carleo2014 ; krutitsky2014 , as can be realized in ultracold-atom experiments greiner2002b ; cheneau2012 ; langen2013 ; geiger2014 . We start from the ground state for some initial value of the interaction parameter $(U/J)_{0}$ and let the system evolve with a different value of $U/J$. In the following, we consider a variety of quenches, spanning the phase diagram, see arrows on Fig. 1(a). We study the spreading of both phase and density fluctuations, via the connected correlation functions $G_{1}(R,t)=\langle\hat{a}^{\dagger}_{R}(t)\hat{a}_{0}(t)\rangle-\langle\hat{a}^{\dagger}_{R}(0)\hat{a}_{0}(0)\rangle$ and $G_{2}(R,t)=g(R,t)-g(R,0)$ with $g(R,t)=\langle\hat{n}_{R}(t)\hat{n}_{0}(t)\rangle-\langle\hat{n}_{R}(t)\rangle\langle\hat{n}_{0}(t)\rangle$. Both can be measured in experiments using time-of-flight and fluorescence microscopy imaging, respectively cheneau2012 ; trotzky2012 ; langen2013 ; geiger2014 .

All the results presented below are obtained using density-matrix renormalization group simulations within the time-dependent matrix-product state ($t$-MPS) representation schollwock2005 ; schollwock2011 ; dolfi2014 . A careful analysis of the numerical cut-offs (high-filling cut-off and bond dimension) has been systematically performed to certify the convergence of the results in all the considered cases. This is particularly critical for quenches in the SF phase where the numerical requirements are most binding note:SupplMat .

Meanfield regime.—

We first consider the meanfield regime in the SF phase, where the numerical results can be compared to analytic predictions. This regime is characterized by a small Lieb-Liniger parameter, $\gamma\equiv U/2J\bar{n}\ll 1$. Figure 2(a) displays the $t$-MPS result for the $G_{2}$ correlation function versus distance ($R$) and time ($t$) for a quench from $(U/J)_{0}=0.2$ to $U/J=0.1$ and $\bar{n}=5$, i.e. from $\gamma_{0}=0.02$ to $\gamma=0.01$ [see red arrow on Fig. 1(a)]. It clearly shows a spike-like structure, characterized by two different velocities. On the one hand, a series of parallel maxima and minima move along straight lines corresponding to a constant propagation velocity $V_{\textrm{m}}$ (the dashed blue lines show fits to two of these minima). On the other hand, the various local extrema start at different activation times $t^{\star}(R)$. The latter are aligned along a straight line with a different slope (solid green line), corresponding a constant velocity $V_{\textrm{\tiny CE}}$. The latter defines the correlation edge (CE) beyond which the correlations are suppressed. Similar results are obtained for all the other quenches in the meanfield regime, as well as for the $G_{1}$ function note:G1mf .

This twofold structure near the CE is readily interpreted using the quasiparticle picture, which we briefly outline here (for details, see Ref. cevolani2018 ): the $G_{1}$ and $G_{2}$ correlation functions are expanded onto the elementary excitations of the system. In the meanfield regime of the BH model, the latter are Bogoliubov quasiparticles with the quasimomentum $k\in[-\pi,+\pi]$ and the dispersion relation

[TABLE]

where $\varepsilon_{k}=4J\sin^{2}(k/2)$ is that of the free-particle tight-binding model. A correlation between two points at a distance $R$ is seeded when two correlated, counter-propagating quasiparticles emanating from the center reach the two points, see Fig. 1(b). The fastest ones are those with the maximum group velocity, V_{\textrm{g}}^{\star}=\underset{k}{\max}\big{(}\hbar^{-1}\partial E_{k}/\partial k\big{)}. It yields the activation time $t^{\star}(R)=R/2V_{\textrm{g}}^{\star}$ and the CE velocity $V_{\textrm{\tiny CE}}=2V_{\textrm{g}}^{\star}$, consistently with the expected Lieb-Robinson bound lieb1972 ; calabrese2006 . More precisely, the correlation at a distance $R$ and a time $t$ is built as a coherent superposition of the contributions of the various quasiparticles. In the vicinity of the CE, only the fastest quasiparticles, i.e. those with a quasimomentum $k$ close to $k^{\star}$, contribute. It creates a sine-like signal at the driving spatial frequency $k^{\star}$, whose extrema move at twice the phase velocity $V_{\varphi}(k)=\hbar^{-1}E_{k}/k$ with $k=k^{\star}$, i.e. $V_{\textrm{m}}=2V_{\varphi}^{\star}$ cevolani2018 . The dispersion around $k^{\star}$ then modulates the sine-like signal by an envelope moving at the CE velocity $V_{\textrm{\tiny CE}}$, see Fig. 1(c). This behavior is reminiscent of the propagation of a coherent wave packet in a dispersive medium brillouin1960 ; lighthill1965 ; born1999 .

To test this picture quantitatively, we have extracted the velocities $V_{\textrm{m}}$ and $V_{\textrm{\tiny CE}}$ from the $t$-MPS results for $G_{2}(R,t)$ by tracking, respectively, the local extrema and the activation time. The results, displayed on Fig. 2(b), show excellent agreement with the theory, i.e. $V_{\textrm{\tiny CE}}\simeq 2V_{\textrm{g}}^{\star}$ and $V_{\textrm{m}}\simeq 2V_{\varphi}^{\star}$ within the fitting errorbars. This cross-validates the $t$-MPS results in the most-demanding SF, meanfield regime on the one hand and the quasiparticle picture above on the other hand. Note that the $t$-MPS results are numerically exact and include effects beyond the Bogoliubov approximation, such as quasiparticle collisions.

Strongly correlated regime at unit filling.—

We now turn to the strongly correlated regime $\gamma\sim 1$, where the correlation functions cannot be systematically computed. We first scan the after-quench interaction parameter $U/J$ from the SF to the MI, along the Mott-$U$ transition at unit filling [${\overline{n}}=1$, see magenta arrows on Fig. 1(a)]. Note that each quench is performed in a unique phase: for $U/J<u_{\textrm{c}}\simeq 3.3$ (SF regime), we use the initial interaction strength $(U/J)_{0}=1$ while for $U/J>u_{\textrm{c}}$ (MI regime), we start from $(U/J)_{0}=\infty$. Figure 3 shows typical results for the spreading of the $G_{1}$ (upper row) and $G_{2}$ (lower row) correlations for quenches to the SF regime [$U/J=0.5$, Fig. 3(a)], and to the MI regime, both slightly beyond the transition [$U/J=8$, Fig. 3(b)], and deep in the MI regime [$U/J=24$, Fig. 3(c)]. In all cases, at the notable exception of $G_{2}$ deep in the MI phase [Fig. 3(c2), see discussion below], we find a twofold spike-like structure. The velocities $V_{\textrm{m}}$ and $V_{\textrm{\tiny CE}}$, extracted as before, are plotted on Fig. 3(d), showing similar results for $G_{1}$ and $G_{2}$. This is consistent with the prediction that these velocities are characterized by the spectrum, irrespective of the observable cevolani2018 .

In the SF regime, $U/J<u_{\textrm{c}}$, the results compare very well with the predictions $2V_{\varphi}^{\star}$ and $2V_{\textrm{g}}^{\star}$ as found from the Bogoliubov dispersion relation (2) [see, respectively, the dashed blue and solid green lines on Figs. 3(d1) and (d2)]. Quite surprizing, the agreement is fair up to the critical point where $\gamma_{\textrm{c}}\simeq 1.6$, far beyond the validity condition of the Bogoliubov theory ($\gamma\ll 1$). In fact, when $U/J$ increases from the meanfield regime, the momentum $k^{\star}$ decreases down to the phonon regime, $k\ll\pi$, and the precise $k$-dependence of the dispersion relation beyond this regime becomes irrelevant. Moreover, the physics being dominated by long wavelength excitations, the lattice discretization in Eq. (1) may be disregarded and the BH model maps onto the continuous Lieb-Liniger model note:SupplMat . The latter is integrable by Bethe ansatz (BA) lieb1963a ; lieb1963b . It yields the sound velocity $V_{\textrm{s}}\simeq 2{\overline{n}}\sqrt{\gamma}\left(1-\sqrt{\gamma}/4\pi\right)$, to lowest order in the weak-$\gamma$ expansion. Up to the critical point, the beyond-meanfield correction, $\sqrt{\gamma}/4\pi$, is less than $10\%$, which explains the good agreement between the numerics and the analytic formula. At the critical point, the numerical results for $V_{\textrm{m}}$ and $V_{\textrm{\tiny CE}}$ are consistent with the exact BA value $2V_{\textrm{s}}\simeq 4.6$ note:Vscrit .

The spreading velocities $V_{\textrm{m}}$ and $V_{\textrm{\tiny CE}}$ are continuous at the Mott-$U$ transition, and do not show any critical behavior. Right beyond the critical point, they are still nearly equal and we can hardly distinguish two features from the numerics up to $U/J\simeq 6$. Deeper in the MI phase, however, we recover two distinct features and two different velocities. Contrary to the SF regime, here we find $V_{\textrm{m}}>V_{\textrm{\tiny CE}}$. These results are readily interpreted from the quasiparticle picture. Deep enough in the MI phase, $U/J\gtrsim 6$, the low-energy excitations are doublon-holon pairs, characterized by the dispersion relation barmettler2012 ; Ejima2012

[TABLE]

The comparison between the spreading velocities $V_{\textrm{m}}$ and $V_{\textrm{\tiny CE}}$ fitted from the $t$-MPS results and the characteristic values $2V_{\varphi}^{\star}$ and $2V_{\textrm{g}}^{\star}$, found from Eq. (3), yields a very good agreement, within less than $5\%$ for $G_{1}$ and $9\%$ for $G_{2}$ [see Figs. 3(d1) and (d2) respectively]. The quantitative agreement between the $t$-MPS results and the theoretical predictions for the $G_{1}$ correlations persists up to arbitrary values of $U/J$. This validates the quasiparticle analysis also in the strong-coupling regime.

Yet, the $G_{2}$ correlations behave differently. For intermediate interactions, $6\lesssim U/J\lesssim 9$, we find a twofold structure consistent with that found for $G_{1}$. The signal for $G_{2}$ blurs when entering deeper in the MI regime, and we are not able to identify two distinct features for $U/J\gtrsim 9$. To understand this behavior, one may resort on a strong-coupling ($U\gg J$) expansion of the correlation functions. In contrast to $G_{1}$, the $G_{2}$ function cannot be cast into the generic form analyzed in Ref. cevolani2018 . Instead, combining Jordan-Wigner fermionization and Fermi-Bogoliubov theory barmettler2012 ; note:SupplMat , one finds $G_{2}(R,t)\simeq-2|g_{2}(R,t)|^{2}$ with

[TABLE]

For $U\gg 2(2\bar{n}+1)J$, the doublon-holon pair dispersion relation (3) reduces to $2E_{k}\simeq U-2(2\bar{n}+1)J\cos(k)$. Owing to the square modulus in the formula $G_{2}(R,t)\simeq-|g_{2}(R,t)|^{2}$, we immediately find that the Mott gap $U$ becomes irrelevant and we are left with the effective dispersion relation $2\tilde{E}_{k}\simeq-2(2\bar{n}+1)J\cos(k)$. On the one hand, the group velocity is not affected and we find the maximum value $2V_{\textrm{g}}^{\star}\simeq 2(2\bar{n}+1)J/\hbar$ at $k^{\star}\simeq\pi/2$. The value $2V_{\textrm{g}}^{\star}=6J/\hbar$ found for ${\overline{n}}=1$ is in excellent agreement with the value of $V_{\textrm{\tiny CE}}$ fitted from the $G_{2}$ function deep in the MI phase, see Fig. 3(d2). On the other hand, the corresponding effective phase velocity vanishes, $2\tilde{V}_{\varphi}^{\star}\simeq 0$. This is consistent with the disappearance of the spike-like structure observed in the $t$-MPS calculations for $G_{2}$ deep in the MI phase note:ZeroPhaseVelocity . In addition, the first-order correction to the leading strong-coupling term, relevant for moderate values of $U/J$, sustains a double structure with $V_{\textrm{g}}^{\star}\neq V_{\varphi}^{\star}$. The latter is consistent with the observation of two distinct spreading velocities, $V_{\textrm{\tiny CE}}\neq V_{\textrm{m}}$, closer to the Mott-$U$ transition note:SupplMat .

Strongly interacting superfluid regime.—

We finally consider the strongly interacting regime of the SF phase, corresponding to $\gamma\gg 1$ and ${\overline{n}}\notin\mathbb{N}$. In this regime, the Tomonaga-Luttinger liquid (TLL) theory accurately describes the low-energy physics of the BH model at equilibrium, including the Mott-$\delta$ transition, see for instance Refs. cazalilla2011 ; haller2010 ; boeris2016 . The TLL theory considers an effective harmonic fluid, characterized by a single characteristic velocity, namely the sound velocity $V_{\textrm{s}}$.

In contrast, our $t$-MPS simulations in the strongly interacting SF regime clearly show beyond TLL physics. We have computed the spreading of correlations for a large value of the after-quench interaction parameter, $U/J=50$, and varying the filling ${\overline{n}}$ up to the Mott-$\delta$ transition at ${\overline{n}}=1$ [see pink arrow on Fig. 1(a)]. The spreading velocities $V_{\textrm{\tiny CE}}$ (green diamonds) and $V_{\textrm{m}}$ (blue disks), found from fits to the two-body correlation function $G_{2}(R,t)$, are shown on Fig. 4. They show clear deviations from twice the sound velocity of the BH model in the strongly interacting limit, 2V_{\textrm{s}}\simeq(4J/\hbar)\sin(\pi{\overline{n}})\big{[}1-(8J/U)\cos(\pi{\overline{n}})\big{]} (orange dotted line and squares) note:VsSISF . Moreover, the emergence of two different characteristic velocities, $V_{\textrm{\tiny CE}}\neq V_{\textrm{m}}$, indicates that the TLL approach is insufficient to describe the spreading of correlations, even upon renormalization of the effective TLL parameters. Note that the two velocities become nearly equal in the vicinity of the Mott-$\delta$ transition and reach the value $V_{\textrm{\tiny CE}}\simeq V_{\textrm{m}}\simeq 6J/\hbar$. This is consistent with the disappearance of the twofold structure and the value found for $V_{\textrm{\tiny CE}}$ deep in the MI phase at ${\overline{n}}=1$, see Fig. 3(d).

Conclusions.—

In summary, working within the case study of the Bose-Hubbard chain and using a numerically-exact many-body approach, we have presented evidence of a universal twofold dynamics in the spreading of correlations. The latter is characterized by two distinct velocities, corresponding to the spreading of local maxima on the one hand and to the CE on the other hand. This has been found in all the phases of the model. Exceptions appear only in a few cases, for instance (i) for specific observables in specific regimes, or (ii) when the two velocities happen to be equal, as found at the Mott critical points for instance.

Our predictions are directly relevant to quench experiments on ultracold Bose gases in optical lattices, where the dynamics of one-body and two-body correlation functions can be observed on space and time scales comparable to our simulations lewenstein2007 ; bloch2008 ; cheneau2012 ; trotzky2012 ; NaturePhysicsInsight2012bloch . Importantly, while in most experiments and numerics the CE is infered from the behavior of the correlation maxima, our results show that the two must be distinguished. This is expected to be a general feature of short-range systems and should be relevant to models other than the sole BH model.

Moreover, our study may be extended to long-range systems, such as spin models as realized in trapped-ion experiments jurcevic2014 ; richerme2014 . While the notions of a maximum group velocity and phase velocity may break down in such systems, the meanfield theory also predicts a twofold dynamics cevolani2018 . In this case, it is characterized by the coexistence of super-ballistic and sub-ballistic signals. The results of the present paper suggest that the twofold structure of the correlation function may survive in strongly correlated regimes also for long-range systems. The demonstration of this effect would shed light on the still debated scaling of the light cone in long-range systems.

This research was supported by the European Commission FET-Proactive QUIC (H2020 grant No. 641122). The numerical calculations were performed using HPC resources from CPHT and GENCI-CCRT/CINES (Grant No. c2017056853), and make use of the ALPS library dolfi2014 . We are grateful to the CPHT computer team for valuable support.

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) schollwock2005 ; schollwock2011 ; dolfi2014 . 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 wolf2008 ; eisert2010 , 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 out 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}\rangle$, where $n_{R}\in\mathbb{N}$, which is infinite. However, the probability distribution of the lattice-site 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_{\textrm{\tiny max}}$. It is important to note that, in some cases, the value of $n_{\textrm{\tiny max}}$ needs to be significantly much larger than the average filling ${\overline{n}}$ and its fluctuations. This observation is consistent with analyses of truncated Bose-Hubbard models in quantum Monte Carlo simulations kashurnikov1998zero .

The SF meanfield regime, which corresponds to a high filling factor $\bar{n}$ and the gapless dispersion relation, has the most binding criteria. We found that a good estimator for $n_{\textrm{\tiny max}}$ is given by the condition $1-\sum_{n=0}^{n_{\textrm{\tiny max}}}P(n)\lesssim 10^{-2}$, where $P(n)$ is the probability that $n$ bosons occupy a given lattice site. In the SF meanfield regime, the probability distribution is nearly Poissonian, $P(n)\simeq\bar{n}^{n}e^{-\bar{n}}/n!$. For instance, for the filling factor ${\overline{n}}=5$ used for the data of Fig. 2, it yields $n_{\textrm{\tiny max}}\gtrsim 12$. For the strongly correlated SF regime at ${\overline{n}}=1$ considered for Fig. 3(a), the density fluctuations are significantly suppressed and using the same condition as previously leads to $n_{\textrm{\tiny max}}=5$. For the MI phase at $\bar{n}=1$ and moderate values of $U/J$ ($15\geq U/J\geq u_{\mathrm{c}}$) considered for Fig. 3(b), we kept $n_{\textrm{\tiny max}}=5$. Deep in the MI phase ($U/J\geq 15$), truncating the local Hilbert space to $n_{\textrm{\tiny 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 factor $\bar{n}<1$ and the large value of the interaction parameter $U/J$, the above condition also yields $n_{\textrm{\tiny max}}=2$, as used for Fig. 4. In all cases, we have checked that the numerics are converged for these values of $n_{\textrm{\tiny max}}$.

Bond dimension.—

Within the MPS approach, the many-body state for a $M$-site lattice is represented in the tensor network form

[TABLE]

where $n_{j}$ spans a local Hilbert space basis. For the BH model, it corresponds to a Fock basis truncated at $n_{\textrm{\tiny max}}$. For each value of $n_{j}$, the quantity $A^{n_{j}}[j]$ is a $\chi_{j-1}\times\chi_{j}$ matrix, where $\chi_{j}$ is the rank associated to the Schmidt matrix when applying the $j$-th singular value decomposition schollwock2011 . The bond dimension $\chi$ is defined as the maximum rank, $\chi=\mathrm{max}_{j}\left(\chi_{j}\right),\leavevmode\nobreak\ j\in[0\ldots M]$. Note that for open-boundary conditions, the quantities $A^{n_{1}}[1]$ and $A^{n_{M}}[M]$ are actually a row vector and a column vector, respectively, i.e. $\chi_{0}=\chi_{M}=1$.

In the numerics, the maximum value of $\chi$ is chosen sufficiently large so that the truncation does not affect the results. In practice, the calculations are run for several values of $\chi$ up to convergence of the correlation function $G_{1}(R,t)$ or $G_{2}(R,t)$. The required value of $\chi$ 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 meanfield regime [Figs. 2(a) and S1], we used the values $\chi=300$ and $\chi=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 ${\overline{n}}=1$ [Fig. 3(a)], we used $\chi=300$ for both correlation functions. A similar value of $\chi$ was considered for moderate values of $U/J$ in the MI phase at $\bar{n}=1$ [Fig. 3(b)]. Deep in the MI phase [Fig. 3(b)], the bond dimension can be significantly decreased and we consider $\chi=100$. Finally, in the SF strongly interacting regime at $U/J=50$, we found that the value $\chi=100$ is enough.

S2 One-body correlation function $\mathbf{G_{1}(R,t)}$ in the meanfield regime

In the analysis of the SF meanfield 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_{\textrm{\tiny CE}}$) and maxima ($V_{\textrm{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{\overline{n}}/J)_{0}=1$ to $U{\overline{n}}/J=0.5$, and $\bar{n}=5$. The fits to the correlation edge and to the maxima yield the velocities $V_{\textrm{\tiny CE}}=(4.4\pm 0.3)\,J/\hbar$ and $V_{\textrm{m}}=(3.3\pm 0.2)\,J/\hbar$, 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 cevolani2018 . 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 bravyi2006 .

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,

[TABLE]

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=\hbar^{2}/2ma^{2}$ and $U=\hbar^{2}c/ma$. The density of the LL model is $\rho\equiv N/L={\overline{n}}/a$, where ${\overline{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 lieb1963a ; lieb1963b . All the thermodynamic quantities at zero temperature can be written as universal functions of the Lieb-Liniger parameter $\gamma=c/\rho$ and the dimensionless quantity $e(\gamma)=E_{0}/Nn^{2}$, where $E_{0}$ is the ground state energy. For instance, the macroscopic sound velocity lieb1963b reads as

[TABLE]

Using the small $\gamma$ expansion, $e(\gamma)=\gamma\left[1-(4/3\pi)\sqrt{\gamma}\right]$, one then finds

[TABLE]

valid in the weakly-interacting regime, $\gamma\ll 1$. Finally, using the correspondance between the parameters of the BH and LL models, one finds

[TABLE]

and $\gamma=U/2J{\overline{n}}$.

S4 Two-body correlation function $\mathbf{G_{2}(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\gg J$ and ${\overline{n}}=1$), we compute the function $G_{2}(R,t)$, working along the lines of Ref. barmettler2012 . 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

[TABLE]

with

[TABLE]

and the excitation spectrum is $2E_{k}\simeq\sqrt{\left[U-2J(2\bar{n}\!+\!1)\cos(k)\right]^{2}+16J^{2}\bar{n}(\bar{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\gg 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)\simeq-2|g_{2}(R,t)|^{2}$,

[TABLE]

Moreover, the excitation spectrum may be expanded in powers of $J/U$. Up to first-order, it yields $2E_{k}\simeq U-2J(2\bar{n}+1)\cos(k)$. The gap term $\textrm{e}^{iUt}$ 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\tilde{E}_{k}=-2J(2\bar{n}+1)\cos(k)$, we then find $G_{2}\simeq-2|g_{2}(R,t)|^{2}$ with

[TABLE]

The integral may be evaluated using the stationary phase approximation. In the infinite time and distance limit along the line $R/t=\mathrm{cst}$, the integral in Eq. (S10) is dominated by the momentum contributions with a stationary phase (sp), i.e. $\partial_{k}(2\tilde{E}_{k}t\pm kR)=0$ or, equivalently, $2\tilde{V}_{\mathrm{g}}(k_{\mathrm{sp}})=\pm R/t$ where $\tilde{V}_{\mathrm{g}}=\partial_{k}\tilde{E}_{k}$ is the group velocity of the effective excitation spectrum. Since the latter is upper bounded by the value $\tilde{V}_{\mathrm{g}}^{*}=\mathrm{max}(\tilde{V}_{\mathrm{g}})=J(2\bar{n}+1)$, it has a solution only for $R/t<2\tilde{V}_{\mathrm{g}}^{*}$. We then find

[TABLE]

with $\sigma=\mathrm{sgn}\left(\partial_{k}^{2}\tilde{E}_{k_{\mathrm{sp}}}\right)$. For both the real and imaginary parts of $g_{2}(R,t)$, the correlations are activated ballistically at the time $t=R/2\tilde{V}_{\mathrm{g}}^{*}$. It defines a linear correlation edge (CE) with velocity $V_{\textrm{\tiny CE}}=2\tilde{V}_{\mathrm{g}}^{*}$. In addition, Eq. (S11) also yields a series of local maxima, defined by the equation $2\tilde{E}_{k_{\mathrm{sp}}}t-k_{\mathrm{sp}}R=\mathrm{cst}$. In the vicinity of the CE cone, these maxima (m) propagate at the velocity $V_{\textrm{m}}=2\tilde{V}_{\varphi}^{*}=2\tilde{E}_{k^{*}}/k^{*}$, i.e. twice the phase velocity at the maximum of the group velocity, $k^{*}$.

Hence, the real and imaginary parts of $g_{2}(R,t)$ both display a twofold structure with a CE velocity $2\tilde{V}_{\mathrm{g}}^{*}=2J(2\bar{n}+1)$ and a velocity of the maxima $2\tilde{V}_{\varphi}^{*}=0$, as shown on Figs. S2(a) and (b). In contrast, $G_{2}(R,t)$, does not display the twofold structure. This is because it is the sum of the squares of the two latter contributions [see Eq. (S11)], which are shifted by half a period and cancel each other. It thus gives a single cone structure, characterized by the sole CE velocity $2\tilde{V}_{\mathrm{g}}^{*}$, as shown on Fig. S2(c).

Quench into the Mott insulator phase for moderate $U/J$.—

For moderate values of $U/J$, still in the MI phase, the second term in the right-hand-side of Eq. (S6), $|\bar{g}_{2}(R,t)|^{2}$, becomes relevant. Using again the stationary-phase approximation for $\bar{g}_{2}(R,t)$, we find

[TABLE]

with $\sigma^{\prime}=\mathrm{sgn}\left(\partial_{k}^{2}E_{k_{\mathrm{sp}}}\right)$ and $E_{k}$ the excitation spectrum given at Eq. (3). Using the same argument as above, we find that $\bar{g}_{2}(R,t)$ shows a twofold structure characterized by, now, the CE velocity $2V_{\textrm{g}}^{*}=2\mathrm{max}\left(\partial_{k}E_{k}\right)$ but the velocity of the maxima $2V_{\varphi}^{*}=2E_{k^{*}}/k^{*}\neq 0$. Since there is a single contribution here, the quantity $|\bar{g}_{2}(R,t)|^{2}$ displays a twofold structure with the same characteristic velocities. More precisely, both the length and time scales of the oscillations are divided by two but the velocities are not affected.

For a quench into the MI phase at a moderate value of $U/J$, both $|g_{2}(R,t)|^{2}$ and $|\bar{g}_{2}(R,t)|^{2}$ contribute to the two-body correlation function $G_{2}(R,t)$. While the $|g_{2}(R,t)|^{2}$ contribution is characterized by the sole CE velocity $2V_{\textrm{g}}^{*}$, the $|\bar{g}_{2}(R,t)|^{2}$ contribution provides the double structure observed on $G_{2}$ for $6<U/J<10$ in the $t$-MPS calculations.