Hybrid semiclassical theory of quantum quenches in one dimensional systems
HHybrid semiclassical theory of quantum quenches in one dimensional systems
C˘at˘alin Pa¸scu Moca,
1, 2
M´arton Kormos, and Gergely Zar´and BME-MTA Exotic Quantum Phase Group, Institute of Physics,Budapest University of Technology and Economics, H-1111 Budapest, Hungary Department of Physics, University of Oradea, 410087, Oradea, Romania BME-MTA Statistical Field Theory Research Group, Institute of Physics,Budapest University of Technology and Economics, H-1111 Budapest, Hungary (Dated: March 20, 2019)We develop a hybrid semiclassical method to study the time evolution of one dimensional quantumsystems in and out of equilibrium. Our method handles internal degrees of freedom completelyquantum mechanically by a modified time evolving block decimation method, while treating orbitalquasiparticle motion classically. We can follow dynamics up to timescales well beyond the reach ofstandard numerical methods to observe the crossover between pre-equilibrated and locally phaseequilibrated states. As an application, we investigate the quench dynamics and phase fluctuations ofa pair of tunnel coupled one dimensional Bose condensates. We demonstrate the emergence of soliton-collision induced phase propagation, soliton-entropy production and multistep thermalization. Ourmethod can be applied to a wide range of gapped one-dimensional systems.
Fundamental questions concerning the coherent timeevolution, relaxation and thermalization of isolated quan-tum systems have been brought into the focus of atten-tion by recent progress in experimental techniques [1–4]. Experiments on cold atomic gases allow us to engi-neer a broad range of lattice and continuum Hamiltoni-ans in a controlled fashion, and to monitor the coher-ent time evolution of these systems through measuringmulti-point correlation functions [5], accessing the quan-tum state via site-resolved quantum microscopy [6], andeven measuring the entanglement properties of the sys-tem [7]. These experiments as well as ongoing matterwave interferometry [8] experiments call for the devel-opment of new analytical and numerical methods thatare able to describe non-equilibrium dynamics in closedinteracting quantum systems, and address fundamentalquestions of non-equilibrium thermodynamics, such asthermalization, entropy production, or the fate of pre-thermalized states.Here we propose a novel method we dub “semi-semiclassical” (sSC), able to reach time scales much be-yond conventional methods [9–11] and to capture thefundamental phenomena of pre-thermalization [12, 13] aswell as local equilibration in great detail. Our method hy-bridizes a semiclassical (SC) approach with time-evolvingblock decimation (TEBD): we compute the time evolu-tion of the internal degrees of freedom completely quan-tum mechanically using TEBD, while treating the orbitalmotion semiclassically. As a proof of principle, we use ourmethod to describe quantum quenches in the sine-Gordonmodel, relevant for two coupled 1D quasicondensates (seeFig. 1.b), and studied intensely both theoretically [14–18] and in matter wave interference experiments usingnanofabricated atom chips [3, 19]. Our sSC method isdemonstrated to capture multistep thermalization andsoliton-collision induced entropy production efficiently,and is applicable to one dimensional gapped systems hav- ing stable quasiparticle excitations [20].Our method is based on a semiclassical (SC) approach,originally developed to study the dynamics of finite tem-perature systems in equilibrium [21–23] and later suc-cessfully applied to quantum quenches [24–27], i.e. tonon-equilibrium situations in which the system evolvesunitarily starting from some prepared initial state. TheSC method is applicable in gapped systems wheneverthe quasiparticles’ Compton wavelength (or the thermalwavelength) is much shorter than their average separa-tion, d = ρ − , with ρ the quasiparticle density. If the en-ergy of the initial state is not too high, it acts as a weaksource of almost pointlike quasiparticles, following clas- t i m e space(a) (b)(c) σ = + σ = − σ = − ϕ ( x ) = ϕ ( x ) − ϕ ( x ) space ++ − ϕ ( x ) ϕ ( x ) x FIG. 1. (Color online) (a) Space-time diagram of quasi-particle trajectories. Labels σ = ± represent internal soli-ton/antisoliton quantum numbers. Collisions are describedby the quantum mechanical S -matrix (solid black dot). (b)One dimensional Bose condensates coupled by quantum tun-neling. Dynamics of the relative phase is described by thesine–Gordon model. (c) The quench creates soliton-antisolitonpairs at t = 0. Arrows indicate directions of propagation. a r X i v : . [ c ond - m a t . s t a t - m ec h ] M a r sical trajectories [28, 29] (see Fig. 1.a). In the standardSC treatment, quasiparticles are furthermore assumed tomove slowly and thus collisions are described by a univer-sal purely reflective scattering matrix [30]. This univer-sal SC (uSC) approach permits the derivation of preciousanalytical results [22, 23, 26, 27], but also suffers from ar-tifacts [23, 27]; certain correlation functions and expecta-tion values fail to decay and internal degrees of freedomremain just locally entangled. Our method eliminates allthese shortcomings [31]. Model.—
The sine-Gordon model is defined as H = (cid:126) c (cid:90) d x (cid:20) πK Π( x ) + K π [ ∂ x ϕ ( x )] (cid:21) − (cid:90) d x cos[ ϕ ( x )] , (1)with Π( x ) denoting the field conjugate to ϕ ( x ), and c the speed of sound. In case of coupled condensates, ϕ ( x ) = ϕ ( x ) − ϕ ( x ) represents the relative phase ofthe condensates, K stands for the Luttinger parameterof each condensate, and the cosine term accounts forJosephson tunneling between them [32]. For locally in-teracting bosons K ≥
1, with K (cid:29) K > / , and for K < / σ = ± corresponding to changes ϕ → ϕ ± π , re-spectively. For K > / f ( v ). The precise form of f ( v ) depends ondetails of the quench protocol, but turns out to be unim-portant in the present calculation. The initial state con-sidered here resembles to the ones appearing in previousstudies of split condensates [34, 35] as well as other sys-tems [16, 36–41], where the initial state was taken to bea coherent superposition of uncorrelated zero momentumquasiparticle pairs. Within our sSC approach, however,phase coherence of the pairs does not play any role. The semi-semiclassical (sSC) method.—
In the SCapproach, the quantum mechanical average is replaced byan ensemble average over initial positions and velocitiesof solition-antisoliton pairs. The corresponding semiclas-sical configurations consist of space-time diagrams withpairs of straight lines (kink trajectories) originating fromthe t = 0 axis, distributed independently from each other,uniformly in space with random slopes corresponding tothe distribution f ( v ) (see Fig. 1.a), and an initial distri-bution of paired charges.Since quantum mechanical effects become relevant onlyat collisions, we can approximately factorize the manybody wave function into an orbital and a charge part ofthe form | Ψ( x , σ , t ) (cid:105) ≈ S| Ψ orb ( x , t ) (cid:105) ⊗ | χ ( σ , t ) (cid:105) , where | Ψ orb ( x , t ) (cid:105) = (cid:90) d N x N (cid:89) j =1 δ ( x j − x j − v j t ) | x , . . . x N (cid:105) , | χ ( σ , t ) (cid:105) = (cid:88) σ ,σ ,... A σ σ ...σ N | σ σ . . . σ N (cid:105) . (2)Here x j denotes the initial position of the j th kink ofvelocity v j and topological charge σ j , and S stands forsymmetrization.In the original uSC approach charges are treated classi-cally: quasiparticles scatter as impenetrable billiard ballswhile preserving their charges. This follows from the per-fectly reflective universal low energy two-body S -matrix, S ˜ σ ˜ σ σ σ = ( − δ σ ˜ σ δ σ ˜ σ (notice the labeling convention),describing the scattering of quasiparticles with vanishingmomenta. Using this perfectly reflective S-matrix allowsone to obtain a number of exact results for thermal gasesand near adiabatic quantum quenches [22, 23, 26, 27].However, as soon as the quench is not slow enough, fasterquasiparticles are inevitably created. Collisions involv-ing these faster particles are not captured by the uni-versal scattering matrix, and the true velocity depen-dent scattering matrix must be used (cf. Fig. 1.a). Thesine–Gordon model is integrable and this two-particle S -matrix is exactly known [42, 43]. The matrix elements S + − + − = S R ( v , v ) and S − ++ − = S T ( v , v ), in particular,describe reflection and transmission, respectively. Theysatisfy | S T | + | S R | = 1 , and for small velocities trans-mission vanishes as | S T ( v , v ) | ∝ ( v − v ) . We can incorporate the non-trivial S -matrix in two dif-ferent ways. In the first method, charges are still treatedclassically, but at each collision either a perfect transmis-sion or perfect reflection takes place with probabilitiesgiven by the modulus square of the S -matrix elements.This method neglects interference effects but can be im-plemented as a classical Monte–Carlo simulation. A morerefined approach is to treat the charge part of the wavefunction in Eq. (2) in a fully quantum mechanical man-ner with an MPS-based method. As the quench protocolgenerates neutral kink pairs with equal probability am-plitude, the initial wave function is a dimerized state, i.e.a product of Bell pairs: | χ ( σ , t = 0) (cid:105) = N/ (cid:89) j =1 | + (cid:105) j − |−(cid:105) j + |−(cid:105) j − | + (cid:105) j √ . (3)The charge wave function evolves only through collisionsand it is frozen between collisions. If line j and j + 1intersect at time t k , the change of the wave function is | χ ( σ , t k, + ) (cid:105) = ˆ S j,j +1 ( t k ) | χ ( σ , t k, − ) (cid:105) , (4)where ˆ S j,j +1 ( t k ) acts nontrivially only on charges j and j + 1 via the 2-body S -matrix evaluated at the velocitiesof the colliding quasiparticles. In this way we mapped thedynamics of the charges to that of an effective quantumspin chain. The time evolution operator of this spin chaindepends on the underlying semiclassical trajectories andis given as a product of local unitary two-body operators,efficiently treated by an MPS-based algorithm. Phase distribution.—
The space resolved relativephase ϕ ( x ) of the condensates can be measured directlyvia interferometry experiments [3, 5, 19, 44]. We firstuse the classical Monte Carlo approach to determinethe full phase distribution after the quench, a quantityalso analyzed experimentally [5]. in the experiment ofRef. [5]. Semiclassically, each kink is a domain wall thatseparates two domains with a phase difference of ± π. Domains separated by s kinks have a phase difference∆ ϕ = 2 π (cid:80) si =1 σ i . The phase in a given domain beingconstant, its distribution function is a sum of weightedDirac-delta peaks located at integer multiples of 2 π , P ( ϕ, t ) = (cid:88) n ∈ Z c πn ( t ) δ (cid:0) ϕ − πn (cid:1) . (5)In experiments, these delta peaks get broadened by quan-tum and thermal fluctuations. To take this into account,we estimated the phase fluctuations (cid:104) ϕ ( x ) (cid:105) around theminima of the cosine potential for typical parameters [43],and broadened the delta-functions accordingly.The resulting phase distribution P ( ϕ, t ) is shownin Fig. 2.a at different instances for a typical smallquench. For the momentum distribution of quasiparti-cles, f ( p ), we used the ansatz f ( p ) ∝ p exp( − p /p ),motivated by overlap expressions in the transverse fieldIsing model [45]. Immediately after the quench, there isa single central peak at ϕ = 0 and c (0) = 1, but sidepeaks emerge as the system evolves in time. The twoside peaks at ϕ = ± π emerge already at the uSC level:initially, domains of phase ϕ = ± π grow ballistically be-tween separating soliton-antisoliton pairs (see Fig. 1.c),therefore c decreases while c ± π increase linearly in time.Within the uSC approach, the weights c ± π saturate at c ± π → / c → / pre-thermalized ’ behavior by sSC at in-termediate times in simulations with small quasiparticle (cid:1)(cid:2) (cid:1) (cid:3) (cid:1) (cid:2) (cid:4) (cid:3) (cid:2) (cid:3) (cid:1) (cid:1) (cid:2) (cid:3) (cid:1) (cid:3) (cid:2) (cid:3) (cid:1) (cid:4) (cid:2) (cid:1)(cid:2) (cid:5) (cid:2)(cid:3) (cid:2)(cid:3) (cid:5) (cid:2) (cid:6) (cid:6) (a) (b) FIG. 2. (a) Phase distributions P ( ϕ, t ) at different times fora system of N = 24 kinks with (cid:104)| v |(cid:105) = 0 . c and Luttingerparameter K = 9 .
82. (b) Time dependence of the weights c πn defined in Eq. (5). velocities (cid:104)| v |(cid:105) /c (cid:46) .
05. However, after a fast relaxation,collisions start to dominate beyond the collision time, t (cid:38) τ = ρ (cid:82) ∞ d v vf ( v ) , and phase propagation takesplace due to transmissive collisions, giving rise to do-mains of phases ϕ = ± π, ± π, . . . , and the emergenceof further side peaks (see Fig. 2.b). These are all expectedto vanish as c πn ∼ / √ t for very long times, reflectingthe expected diffusive nature of phase propagation, asindeed supported by our numerical results (see [43]). Entanglement entropy.—
Our method is able to fol-low the propagation and growth of entanglement in thecharge sector. The most widely used entanglement mea-sure is the entanglement entropy of a subsystem A ,defined as the von Neumann entropy of its reduceddensity matrix: S = − Tr A (cid:2) ρ A log ρ A (cid:3) where ρ A =Tr ¯ A | Ψ( x , σ , t ) (cid:105)(cid:104) Ψ( x , σ , t ) | . Here we focus on the simplecase when A is half of the total system. The initial chargewave function | χ (cid:105) has alternating bond entanglement en-tropies of log 2 and 0. At t = 0 + , the spatial extension ofthe entangled bonds is zero, so with probability 1 the cutbetween the two halves of the system falls between entan-gled pairs resulting in zero entanglement entropy. Withinthe uSC approximation, at t >
0, pairs with one memberon the left and another member on the right entangle thetwo halves of the system. Since the uSC S -matrix is fullyreflective, each kink remains maximally entangled withits original partner, but entanglement remains local. Atvery long times, the cut between the two halves of thesystem either falls between pairs or cuts a pair in two withprobability 1/2. Consequently, the entanglement entropysaturates at a value S sat = log 2 / . The full time evolu-tion of the entropy is simple to calculate within the uSCapproach which yields an exponential relaxation to S sat , S uSC ( t ) = (1 − e − t/τ ) ln 2 / , (solid black line in Fig. 3).The saturation of S uSC ( t ) is clearly an artefact of theuSC method. In reality, the topological charges of remotequasiparticles become entangled with time due to trans- t= = S ( t ) = ( l n ( ) = ) hj v ji =c = 0 :
04 (uSC)0.04 (sSC)0.080.11
FIG. 3. Dynamics of charge entanglement entropy after thequench for the same parameters as in Fig. 2. The solid blackline denotes the uSC result with a fully reflective S -matrix. missive scattering processes. The sSC approach capturesthe corresponding entropy production: the entanglemententropy does not saturate after the initial transient, butis found to grow linearly in time without bound in aninfinite system, S ( t ) = α t , with a growth rate α ∝ (cid:104)| v |(cid:105) . Correlation functions and equilibration.—
The sSCmethod is also suitable for computing the out of equi-librium evolution of correlation functions. We considerfirst the expectation values G α ( t ) = (cid:104) e iαϕ ( x,t ) (cid:105) = (cid:104) e iαϕ ( x,t ) e − iαϕ ( x, (cid:105) . The function G ( t ) is essentially thecoherence factor, directly accessible through matter waveinterferometry [3]. The standard uSC approach has beenrecently used to compute G α ( t ) [27], which was foundto decay exponentially to a value cos ( πα ), apart from aprefactor incorporating vacuum fluctuations, now set toone (see Fig. 4.a). This surprising behavior is related tothe fact that in the uSC approach the phase is pinned tothe values ϕ = 0 , ± π . In contrast, within the sSC ap-proach, the phase meanders with time and G α ( t → ∞ )no longer remains finite; after a fast transient describedapproximately by uSC, G α ( t ) is found to decay exponen-tially to zero with a rate depending on (cid:104)| v |(cid:105) /c . The equaltime two-point correlation functions, C α ( x − x (cid:48) ; t ) = (cid:104) e iαϕ ( x,t ) e − iαϕ ( x (cid:48) ,t ) , (cid:105) (6)are also accessible experimentally [46], and the integratedquantity, (cid:82) d x d x (cid:48) C ( x − x (cid:48) ; t ) , is directly related to thecontrast of the interference fringes. Fig. 4.b shows thetime evolution of C α ( x ; t ) for α = 1 /
2. Here again, theuSC approach yields a quick relaxation to a state, wherethe phase is pinned and, accordingly, long-ranged corre-lations persist, C α ( | x | , t → ∞ ) ∼ cos ( πα ).In sharp contrast to traditional uSC, within the sSCwe find that C α ( x − x (cid:48) , ∞ ) = exp( − ( πα ) ρ | x − x (cid:48) | )which formally coincides with the thermal equilibriumcorrelator computed within the standard finite temper-ature SC approach [22]. In our case, however, the finalstate is not thermal: in the spirit of Generalized GibbsEnsemble (GGE) the density ρ is conserved and set by (cid:1) (cid:1) -6 -4 -2 0 2 4 6 (cid:2) (cid:1) (cid:3) (cid:2) (cid:3) (cid:1) (cid:4) (cid:5) (cid:4) (cid:1)(cid:2) (cid:6) (cid:2)(cid:3) (cid:3)(cid:7) (cid:5) (cid:6) (cid:8) (cid:7) (cid:5)(cid:3) (cid:2) (cid:5) (cid:2)(cid:8)(cid:2)(cid:6) (cid:4) (cid:3)(cid:9)(cid:10)(cid:11)(cid:12)(cid:13)(cid:14)(cid:15)(cid:4) (cid:1)(cid:2) (cid:1) (cid:3) (cid:1) (cid:2) (cid:2) (cid:3) (cid:1) (cid:4) (cid:1)(cid:2) (cid:4) (cid:2)(cid:3) (cid:2)(cid:5) (cid:5) (cid:6) (cid:6) (cid:6)(cid:7) (cid:3)(cid:8)(cid:9)(cid:10)(cid:4)(cid:6)(cid:11)(cid:6)(cid:7) (cid:3)(cid:12)(cid:9)(cid:10)(cid:4)(cid:6)(cid:11)(cid:1)(cid:13)(cid:6)(cid:11)(cid:14) (a)(b) FIG. 4. (a) Relaxation of the expectation value (cid:104) e iϕ ( x,t ) / (cid:105) for different average velocities. The solid black line representsthe uSC result (fully reflective S-matrix). (b) Equal time cor-relation function (cid:104) e iϕ ( x,t ) / e − iϕ ( x (cid:48) ,t ) / (cid:105) at different times cal-culated with the sSC method. The black dashed line is thethermal equilibrium result in the uSC approximation. the initial state rather then a single effective temper-ature [47]. The final state we obtain should rather beviewed as a pre-thermalized state, where local phase cor-relations and expectation values of vertex operators lookthermal, but the velocity distribution of the quasiparti-cles remains non-thermal, and the quasiparticles’ averageenergy is not related to their density. The sSC methodthus apparently captures aspects of local phase equilibra-tion and “pre-thermalization”.Equilibration is thus predicted to take place in severalsteps in a split condensate. For small values of (cid:104)| v |(cid:105) /c ,first a quick relaxation occurs to a first pre-equilibratedstate with pinned phase and a non-thermal quasiparticlevelocity distribution, but with kink positions random-ized. Next, the coherent evolution of the charge wavefunction gives rise to phase (quantum) propagation. Atthis stage, phase correlation functions relax to their GGEvalues. To reach a truly thermal state with thermal quasi-particle velocity distribution a further relaxation stepand coupling to some external environment such as thesymmetrical phase mode are ultimately needed in thecoupled condensate experiment. Conclusions.—
The versatile semi-semiclassicalmethod developed here has a broad range of applicabil-ity. It is suitable for studying the dynamics of latticesystems and spin chains [33] as well as continuum 1Dsystems. The differences between these translate into dif-ferences in the nature of quasiparticles and their 2-bodyS-matrix, which makes our method ideal for identifyinguniversal aspects of the dynamics. The semiclassicalaspect allows for heuristic interpretations and relativelylong simulation times, while the quantum description ofthe scattering of quasiparticles allows us to go beyondthe limitations of the standard uSC approach and opensthe way to study the propagation of entanglement.Possible connections with the recent work [48] wouldbe interesting to analyze. Although here we focused onthe non-equilibrium time evolution, the method can alsobe used to investigate dynamical correlation functionsat finite temperature beyond the uSC approximation,and is well-suited to study non-equilibrium dynamics ininhomogeneous 1D systems [49].
Acknowledgements.
We thank J¨org Schmiedmayer,Eugene Demler and Ehud Altman for fruitful discussions.This work was supported by the Hungarian research fundOTKA under grant Nos. K105149 and SNN118028. M.K.was partially supported by a J´anos Bolyai Scholarshipand a Pr´emium Postdoctoral Fellowship of the HAS. [1] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat-tore, Rev. Mod. Phys. , 863 (2011).[2] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature ,900 (2006).[3] S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm,and J. Schmiedmayer, Nature , 324 (2007).[4] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa,B. Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler,and J. Schmiedmayer, Science , 1318 (2012).[5] T. Schweigler, V. Kasper, S. Erne, B. Rauer, T. Langen,T. Gasenzer, J. Berges, and J. Schmiedmayer, ArXive-prints (2015), arXiv:1505.03126 [cond-mat.quant-gas].[6] W. S. Bakr, J. I. Gillen, A. Peng, S. F¨olling, andM. Greiner, Nature , 74 (2009).[7] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin,M. Rispoli, and M. Greiner, Nature , 77 (2015).[8] T. Schumm, S. Hofferberth, L. Andersson, S. Wilder-muth, S. Groth, Bar, J. Schmiedmayer, and P. Kruger,Nat. Phys. , 57 (2005).[9] S. R. White and A. E. Feiguin, Phys. Rev. Lett. ,076401 (2004).[10] U. Schollwock, Annals of Physics , 96 (2011).[11] G. Vidal, Phys. Rev. Lett. , 147902 (2003).[12] J. Berges, S. Bors´anyi, and C. Wetterich, Phys. Rev.Lett. , 142002 (2004).[13] M. Buchhold, M. Heyl, and S. Diehl, Phys. Rev. A ,013601.[14] V. Gritsev, E. Demler, M. Lukin, and A. Polkovnikov,Phys. Rev. Lett. , 200404 (2007).[15] A. Iucci and M. A. Cazalilla, New J. Phys. , 055019(2010).[16] B. Bertini, D. Schuricht, and F. H. L. Essler, J. Stat.Mech. , P10035 (2014).[17] C. Neuenhahn and F. Marquardt, New Journal of Physics , 125007 (2015).[18] G. Delfino and J. Viti, J. Phys. A: Math. Theor. ,084004 (2017).[19] M. Kuhnert, R. Geiger, T. Langen, M. Gring, B. Rauer,T. Kitagawa, E. Demler, D. Adu Smith, and J. Schmied-mayer, Phys. Rev. Lett. , 090405 (2013). [20] The most common methods, TEBD [11] and t-DMRG [9,10], are both limited by finite size effects, and typi-cally work only for short times, preceding thermalization.Quantum Monte Carlo techniques may be a promisingdirection (see G. Carleo et al., (2016) arXiv:1612.06392).[21] S. Sachdev and A. P. Young, Phys. Rev. Lett. , 2220(1997).[22] K. Damle and S. Sachdev, Phys. Rev. Lett. , 187201(2005).[23] A. Rapp and G. Zar´and, Phys. Rev. B , 014433 (2006).[24] H. Rieger and F. Igl´oi, Phys. Rev. B , 165117 (2011).[25] D. Rossini, A. Silva, G. Mussardo, and G. Santoro, Phys.Rev. Lett. , 127204 (2009).[26] S. Evangelisti, J. Stat. Mech. , P04003 (2013).[27] M. Kormos and G. Zar´and, Phys. Rev. E , 062101(2016).[28] P. Calabrese and J. Cardy, Phys. Rev. Lett. , 136801(2006).[29] P. Calabrese and J. Cardy, J. Stat. Mech. , P04010(2005).[30] S. Sachdev and K. Damle, Phys. Rev. Lett. , 943(1997).[31] We follow the nomenclature of Ref. [21]. There low den-sity implies small quasiparticle momenta, p , and since p ∼ (cid:126) , the limits p → (cid:126) → p but require low quasiparticle density.[32] V. Gritsev, A. Polkovnikov, and E. Demler, Phys. Rev.B , 174511 (2007).[33] Detailed iTEBD simulations on the spin S = 1 Heisen-berg chain validity of the sSC approach in the dilute limit[M. Werner, C.P. Moca, M. Kormos, and G. Zar´and, un-published].[34] T. Kitagawa, A. Imambekov, J. Schmiedmayer, andE. Demler, New J. Phys. , 073018 (2011).[35] L. Foini and T. Giamarchi, Phys. Rev. A , 023627(2015).[36] A. Silva, Phys. Rev. Lett. , 120603 (2008).[37] S. Sotiriadis and J. Cardy, Phys. Rev. B , 134305(2010).[38] D. Fioretto and G. Mussardo, New J. Phys. , 55015(2010).[39] J. De Nardis, B. Wouters, M. Brockmann, and J.-S.Caux, Phys. Rev. A , 033601 (2014).[40] B. Pozsgay, J. Stat. Mech. , P06011 (2014).[41] D. Horv´ath, S. Sotiriadis, and G. Tak´acs, Nucl. Phys. B , 508 (2016).[42] A. B. Zamolodchikov and A. B. Zamolodchikov, Annalsof Physics , 253 (1979).[43] C. P. Moca, M. Kormos, and G. Zar´and, Supplementarymaterial .[44] G. B. Jo, J. H. Choi, C. A. Christensen, Y. R. Lee, T. A.Pasquini, W. Ketterle, and D. E. Pritchard, Phys. Rev.Lett. , 240406 (2007).[45] P. Calabrese, F. H. L. Essler, and M. Fagotti, J. Stat.Mech. , P07016 (2012).[46] T. Langen, R. Geiger, M. Kuhnert, B. Rauer, andJ. Schmiedmayer, Nat. Phys. , 1 (2013).[47] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii,Physical Review Letters , 050405 (2007).[48] V. Alba and P. Calabrese, (2016), arXiv:1608.00614.[49] M. Kormos, C.P. Moca, and G. Zar´and, unpublished.[50] In the present calculation the velocity distribution is de-rived from a distribution with relativistic dispersion, but we have checked that the results remain largely unaltered by using different distributions. Supplemental Material
The Hamiltonian and sine–Gordon description of coupled quasicondensates
The Hamiltonian of two coupled 1D condensates is given by [32] H = (cid:88) j =1 , (cid:90) d x (cid:26) (cid:126) m ∂ x ψ † j ( x ) ∂ x ψ j ( x ) + g ψ † j ( x ) ψ † j ( x ) ψ j ( x ) ψ j ( x ) + [ V ( x ) − µ ] ψ † j ( x ) ψ j ( x ) (cid:27) − (cid:126) J (cid:90) d x (cid:104) ψ † ( x ) ψ ( x ) + ψ † ( x ) ψ ( x ) (cid:105) , (S-7)where ψ ( x ) , ψ ( x ) are the bosonic fields in the two condensates, V ( x ) is the longitudinal potential which we neglectin the following, J is the tunnel coupling between the condensates, and g = 2 (cid:126) a s ml ⊥ (cid:18) − . a s l ⊥ (cid:19) − (S-8)is the strength of the atom-atom interaction within each condensate. Here l ⊥ = (cid:112) (cid:126) / ( mω ⊥ ) with ω ⊥ being thefrequency of the radial confinining potential and a s denotes the s -wave scattering length of the atoms. The strengthof interaction is often parameterized by the dimensionless combination γ = mg (cid:126) n . (S-9)In the absence of coupling, the condensates can be described within the bosonization framework by the Hamiltonians H j = (cid:126) c (cid:90) d x (cid:26) πK Π j ( x ) + Kπ [ ∂ x ϕ j ( x )] (cid:27) , (S-10)where [ ϕ j ( x ) , Π j ( x (cid:48) )] = iδ ( x − x (cid:48) ) . The speed of sound, c, and the Luttinger parameter K can be computed from theexact Bethe Ansatz solution of the model. The asymptotic expansions for small and large couplings are K ≈ π √ γ (cid:18) − √ γ π (cid:19) − / ≈ (cid:126) π (cid:114) nmg , c ≈ (cid:114) ngm for γ (cid:46) , (S-11a) K ≈ (1 + 4 /γ ) , c ≈ (cid:126) πn/m for γ (cid:29) . (S-11b)Due to Galilean invariance, cK = (cid:126) nπ/m holds for all γ. The density fluctuations are suppressed at wavelengthssmaller than the healing length which is also used as a short distance cutoff. For small γ it is ξ h = 1 / ( n √ γ ) = (cid:126) / √ mgn ≈ (cid:126) /mc , (S-12)while at strong coupling ξ h ≈ /n. The coupling between the condensates is captured by the term 2∆ cos( ϕ − ϕ ) , where ∆ ≈ (cid:126) Jn for weak coupling, but it can be renormalized at strong interactions. It is convenient to introduce thetotal and relative phase ϕ ± = ϕ ± ϕ and Π ± = (Π ± Π ) / H − = (cid:126) c (cid:90) d x (cid:20) πK Π − ( x ) + K π [ ∂ x ϕ − ( x )] (cid:21) − (cid:90) d x
2∆ cos[ ϕ − ( x )] . (S-13) Fluctuations of the phase
We are interested in the fluctations of ϕ − around the minima of the cosine potential. We will calculate this in theharmonic approximation, i.e. by expanding the cosine up to quadratic order which yields H − ≈ (cid:126) c (cid:90) d x (cid:20) πK Π − ( x ) + K π [ ∂ x ϕ − ( x )] + 2∆ (cid:126) c ϕ − ( x ) (cid:21) , (S-14)a free massive boson theory with mass gap m = (cid:114) (cid:126) πc K = (cid:114) mnc . (S-15)The mode expansion of the fields are ϕ ( x ) = (cid:114) πcK (cid:90) d p π (cid:112) ω ( p ) (cid:104) b p e ipx/ (cid:126) + b † p e − ipx/ (cid:126) (cid:105) , (S-16)Π( x ) = − i (cid:126) (cid:114) Kπc (cid:90) d p π (cid:112) ω ( p ) (cid:104) b p e ipx/ (cid:126) − b † p e − ipx/ (cid:126) (cid:105) (S-17)with ω ( p ) = (cid:112) p c + m c and [ b p , b † p ] = 2 πδ ( p − p (cid:48) ) . The finite temperature equal time correlation function of the phase is (cid:104) ϕ ( x ) ϕ ( x (cid:48) ) (cid:105) T = πcK (cid:90) d p π ω ( p ) e ip ( x − x (cid:48) ) / (cid:126) coth( βω ( p ) / , (S-18)where we used (cid:10) b † p b p (cid:48) (cid:11) = 2 πδ ( p − p (cid:48) ) f T ( p ) with f T ( p ) = 1 / ( e ω ( p ) / ( k B T ) −
1) the thermal Bose–Einstein distributionand β = 1 / ( k B T ) . Quantum fluctuations ( T = 0 ) At zero temperature the integral can be evaluated in closed form, (cid:104) ϕ ( x ) ϕ ( x (cid:48) ) (cid:105) T =0 = πcK (cid:90) d p π ω ( p ) e ip ( x − x (cid:48) ) / (cid:126) = 1 K K [ m c | x − x (cid:48) | / (cid:126) ] = 1 K K (cid:18) ∆ xl ∆ (cid:19) , (S-19)where ∆ x = | x − x (cid:48) | , K ( z ) is the modified Bessel function of the second kind and l ∆ ≡ (cid:126) m c = (cid:114) (cid:126) n m ∆ (S-20)is the Compton wavelength associated with the mass m which physically corresponds to the “healing length of therelative phase”. For small separation (cid:104) ϕ ( x ) ϕ ( x (cid:48) ) (cid:105) T =0 = − K (cid:40) log (cid:18) ∆ x l ∆ (cid:19) (cid:34) (cid:18) ∆ x l ∆ (cid:19) (cid:35) + γ E (cid:41) + O (cid:34)(cid:18) ∆ xl ∆ (cid:19) (cid:35) , (S-21)so it is logarithmically divergent ( γ E is the Euler–Mascheroni constant). By introducing a minimum distance α weobtain (cid:10) ϕ (cid:11) T =0 ≈ − K (cid:20) log (cid:18) α l ∆ (cid:19) + γ E (cid:21) . (S-22)For e.g. l ∆ /α ≈ , (cid:10) ϕ (cid:11) T =0 ≈ . /K. As K > , this means that the quantum fluctuation of the phase satisfies∆ ϕ (cid:46) . π, so the approximation of the cosine by a quadratic potential is justified. t= = c : n -3 -2 -1 c c ' : c ' : c ' : hj v ji =c = 0 : FIG. S-5. (Color online) Time dependence of the weights c πn defined in Eq. (5) in the main text represented on a logarithmicscale. The dotted blue line represents the expected ∼ (cid:112) τ /t dependence corresponding to diffusive behavior in the large timelimit. Thermal fluctuations (low temperature)
Let us now compute the correlation function focusing on the low temperature limit, i.e. when k B T (cid:28) m c . Thethermal contribution is given by the expression (cid:104) ϕ ( x ) ϕ ( x (cid:48) ) (cid:105) therm = πcK (cid:90) d p π ω ( p ) e ip ( x − x (cid:48) ) / (cid:126) (coth( βω ( p ) / − , (S-23)which is perfectly well-behaved for large p. Expanding the hyperbolic cotangent in powers of e − βω ( p ) , changing themomentum integration variable to relativistic rapidity, and shifting the integration contour we arrive at (cid:104) ϕ ( x ) ϕ ( x (cid:48) ) (cid:105) therm = 2 K ∞ (cid:88) n =1 K (cid:16)(cid:112) ( nπq/ K ) + (∆ x/l ∆ ) (cid:17) , (S-24)where q ≡ λ T l ∆ = 2 Kπ m c k B T (S-25)is the ratio of the thermal phase coherence length λ T = 2 (cid:126) n/ ( mk B T ) and l ∆ . For the fluctuations of ϕ this implies (cid:10) ϕ (cid:11) T = (cid:10) ϕ (cid:11) T =0 + (cid:10) ϕ (cid:11) therm ≈ K (cid:34) − log (cid:18) α l ∆ (cid:19) − γ E + 2 ∞ (cid:88) n =1 K (cid:16) n πq K (cid:17)(cid:35) . (S-26)In the low temperature limit πq/ (2 K ) (cid:29) K (cid:16) n πq K (cid:17) = (cid:115) Knq e − nπq/ (2 K ) (cid:32) − K πnq + O (cid:34)(cid:18) Kπnq (cid:19) (cid:35)(cid:33) . (S-27)Each term in the sum in Eq. (S-26) is exponentially suppressed with respect to the previous one, so we can truncatethe series at the first term, leading to (cid:10) ϕ (cid:11) T ≈ K (cid:32) − log (cid:18) α l ∆ (cid:19) − γ E + 2 (cid:115) Kq e − m c / ( k B T ) (cid:33) . (S-28)The calculation is consistent if (cid:10) ϕ (cid:11) T (cid:28) δ -peaks accordingly. For completeness, we also represent in Fig. S-5 the coefficients c πn on a log-logscale which shows a decay ∼ / (cid:112) ( t/τ ) further supporting the diffusive behavior for the propagation of the phase, asdiscussed in the main body of the paper. S-matrix of the sine–Gordon model
The 2-particle S-matrix describing the scattering of two kinks is exactly known [42]. The energy and momentumof the incoming and outgoing particles are conveniently parameterized in terms of the relativistic rapidity as E = mc cosh θ, p = mc sinh θ. The 2-particle S-matrix is given by a four by four matrix in the basis | + + (cid:105) , | + −(cid:105) , | − + (cid:105) , | − −(cid:105) : S = S S T S R S R S T S , (S-29)where due to relativistic invariance all entries depend only on the relative rapidity θ = θ − θ of the two incomingkinks. Here S T is the amplitude of transmission and S R is the amplitude of reflection. The term S accounts for thephase picked up by the wave function upon scattering of two kinks of the same charge: S ( θ ) = − exp (cid:40) − i (cid:90) d tt sinh t ( π − ξ )2 sinh ξt cosh πt sin( θt ) (cid:41) , (S-30)where ξ = π K − . (S-31)The transmission and reflection factors are given by S T ( θ ) = sinh πθξ sinh π ( iπ − θ ) ξ S ( θ ) , (S-32) S R ( θ ) = i sin π ξ sinh π ( iπ − θ ) ξ S ( θ ) . (S-33)and they satisfy | S T | + | S R | = 1 . For small θ they behave as | S T | ∝ θ , | S R | ∝ − θ , thus at small rapidities thescattering of kinks of opposite charges is almost purely reflective. Details of the numerical algorithm
In this section we discuss in more detail the numerical algorithm which is a combination of Monte Carlo (MC)sampling of classical trajectories for the soliton-antisoliton pairs and a propagation in time of the initial wave function | Ψ orb ( x , t ) (cid:105) using the Time-Evolving Block Decimation (TEBD) approach. The process can be divided into three mainsteps as follows:(1) Generation of classical kink configurations – For a given concentration of pairs, we randomly generate theirpositions at t = 0. The kink velocities are generated from a given velocity distribution [50]. Then we constructthe kink configuration (one such configuration is presented in Fig. S-6(a)) which is a space-time diagram in ( x, t )coordinates that displays the classical trajectories of the kinks. We have to keep in mind that at t = 0 two kinks thatform a pair start to move in opposite directions with equal velocities. In this representation the collision of two kinkscorresponds to the intersection of two lines. When constructing the kink configuration we index all the intersectionpoints in the space-time coordinates ( x I , t I ) as well as the corresponding lines. These coordinates are then orderedchronologically. In this way we completely characterize the orbital motion of the kinks, which corresponds to theconstruction or the orbital part of the wave function | Ψ orb ( x , t ) (cid:105) in Eq. (2b). Furthermore (although not displayed inFig. S-6(a)) we impose hard wall boundary conditions such that a kink is perfectly reflected off the walls.(2) Construction of the effective spin model – Quantum effects become relevant only at collision times and they donot affect the orbital motion of the quasiparticles, which allows us to factorize the wave function as in Eq. (2) in themain body of the paper. To be able to propagate the initial wave function we need to map the dynamics of the spinsto an effective spin model. We do that by constructing first zig-zag kink configuration, which consists of relabelingthe lines according to Fig. S-6(b). In this way, the lines are ordered from left to right at any time instance t but, more0 t i m e space t i m e space1 2 3 4 5 6 t i m e space(a) (b) (c) FIG. S-6. (Color online) (a) A typical kink configuration with soliton-antisoliton pairs. Their space-time evolution is describedby pairs of straight lines. Different lines are plotted with different colors. (b) Reindexing the lines after each collision. Thequasiparticle trajectories have now zig-zagged shapes indicated by the colors. With this labeling the collisions always take placebetween neighboring quasiparticles. (c) Propagation of the TEBD solution in the charge sector of the wave function. When twolines intersect in (b), the effective spin chain of the charge sector is acted on by the unitary evolution operator given by thecorresponding S-matrix. importantly, in this representation only neighboring lines intersect. Furthermore, each line j carries a “spin” σ j . Theintersections of classical lines corresponds to a scattering event of the two spins carried by the lines. This is a quantummechanically process that is fully described in terms of the S-matrix, so at this point the model Hamiltonian is notrelevant as all the necessary information that we need is encoded in the two-body S-matrix.(3) Time evolution of the wave function – The mapping of the zig-zag kink configuration to the effective spin chainmodel is displayed in Fig. S-6(c). Here the square regions indicate the interaction between two neighboring spins.Basically, the spin chain is frozen in time in between collisions and the time evolution takes place only when pairs ofspins are scattered at the collision times. This picture allows us to use the TEBD algorithm [10, 11] and propagatethe initial wave function. The t = 0 wave function | χ ( σ , t = 0) (cid:105) given in Eq. (3) is first organised as an MPS statewhich is then evolved in time. The full time evolution operator U ( t ) is the chronological product of unitary two-bodyS-matrix operators that scatter pairs of quasiparticles only at the intersection times t I . Furthermore, the TEBDframework allows us to compute different physical quantities [11]. For example, to compute the time evolution of theentanglement entropy, we first evolve the MPS state up to the desired time t, bipartition the state using the Schmidtdecomposition and then compute the entanglement entropy in the usual way [11].What we have described so far is the evolution of a single kink configuration. To compute the averages discussedin the main body of the paper, we sample in general ∼5