{'fraction_non_alphanumeric': 0.04103927438452259, 'fraction_numerical': 0.013091461507528707, 'mean_word_length': 4.5392896918698185, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum sine-Gordon dynamics in coupled spin chains

The sine-Gordon field theory emerges as the low-energy description in a wealth of quantum many-body systems. Recent efforts have been directed towards realizing quantum simulators of the model, by interfering two weakly coupled one-dimensional cold atomic gases. The weak interactions within the atomic clouds provide a sine-Gordon realization in the semiclassical regime. Furthermore, the complex microscopic dynamics prevents a quantitative understanding of the effective sine-Gordon validity realm. In this work, we focus on a spin ladder realization and observe the emergent sine-Gordon dynamics deep in the quantum regime. We use matrix-product state techniques to numerically characterize the low-energy sector of the system and compare it with the exact field theory predictions. From this comparison, we obtain quantitative boundaries for the validity of the sine-Gordon description. We provide encompassing evidence for the emergent field theory by probing its rich spectrum and by observing the signatures of integrable dynamics in scattering events.


I. INTRODUCTION
In recent years, highly controllable and tunable quantum simulators have been developed on different platforms, including ultracold atoms, trapped ions, and superconducting qubits. These settings enabled the experimental preparation and characterization of stronglycorrelated non-equilibrium states of matter. In general, such states are extremely difficult to characterize by traditional analytical or numerical approaches, which defies a direct verification of the quantum simulation. To overcome this challenge, it is important to identify strongly correlated systems, whose dynamics can also be characterized with conventional theoretical means. In this respect, a prominent example is the sine-Gordon model [1][2][3][4]; a one-dimensional relativistic field theory which describes the low-energy physics of a multitude of experimental systems. The sine-Gordon model is a remarkable example of an integrable field theory [2,5]. Its spectrum features topological excitations, akin to classical solitons, and their bound states, known as breathers. Out of equilbirium, integrability hinders thermalization and many efforts have been made to understand the exotic dynamics of this field theory [6][7][8][9][10][11][12].
The sine-Gordon model captures the low-energy sector of spin chains [13][14][15][16][17][18], spinful cold atom gases [3], specific quantum circuits [19,20] and the interference of two onedimensional quasi-condensates [21,22]. The latter realization of the sine-Gordon model, proposed by Gritsev et al. [21], has been realized with coupled weakly-interacting quasi-condensates on atom chips [23], coupled by a potential barrier of adjustable height. The sine-Gordon Hamiltonian governs the relative phase between the two condensates, which is then probed through matter-wave interferometry [24,25]. Importantly, the mass scale of the theory can be tuned by controlling the potential barrier. For weak interactions within the one dimensional gases, as realized in these experiments, a semiclassical approximation of the sine-Gordon model accurately captures the equilibrium correlation functions of the relative phase [23,26].
Out of equilibrium the situation is less clear, as highenergy excitations beyond the sine-Gordon description are inevitably created. Quench experiments in tunnel coupled quasi-condensates [27] are not described by a semiclassical analysis of the sine-Gordon model [28]. Attempts to include corrections beyond sine-Gordon arising from coupling to the symmetric-phase sector [29] and transverse mode excitations [30] seem not to be sufficient to quantitatively describe the experimental observations [27]. From a theoretical vantage point, both the massless and large mass limits (corresponding to large and weak barrier strengths respectively) become quadratic theories, and their dynamics are readily tractable [31][32][33][34]. Hence, they represent a good starting point for self-consistent Gaussian approximations [35]. Likewise, experiments focusing on these regimes [36][37][38][39][40] are well captured by these approaches. Yet, fascinating physics arises precisely in the challenging regime of intermediate mass scales and far from the semiclassical limit, in which the rich dynamics of the field theory is unveiled. These considerations pose two questions. First, in order to probe strong quantum effects, it will be interesting to focus on the strong-coupling regime beyond the semiclassical limit. Second, assessing the validity regime of the emergent sine-Gordon descriptions in microscopic realizations will be pertinent.
These open questions motivate our work, in which we study the sine-Gordon dynamics emerging from the tunnel coupling between two one-dimensional XXZ spin chains, see Fig. 1, which realizes a lattice version of the original proposal of Ref. [21]. Our motivation to investigate lattice systems is twofold. On the one hand, matrixproduct-state techniques [41][42][43] allow for a precise numerical characterization of the spectrum and the dynamics of the system. Therefore, well-defined boundaries to the effective sine-Gordon dynamics can be obtained. On the other hand, interactions in these settings are strong and result in sine-Gordon realizations deep in the quantum regime, in contrast with experiments on coupled quasi-condensates. A realization of a model closely related to the spin ladder, that consists of coupled Bose-Hubbard chains is realizable with current experimental capabilities [44].
Our discussion is organized as follows: in Section II we introduce the coupled spin chains and review the emergence of the effective sine-Gordon dynamics. Section III provides a short summary of the integrability aspects of the sine-Gordon model, which are then compared in Section IV with the numerical low-energy spectrum of the coupled chains. In Section V we analyze nonequilibrium scenarios and observe key signatures of integrable dynamics in scattering events of wave-packet excitations. The appendices containing more technical aspects then follow the conclusions and outlook VI.

II. THE SINE-GORDON FIELD THEORY ON A SPIN LADDER
The sine-Gordon model is the central character of this work and is governed by the Hamiltonian where the bosonic fields are canonical conjugates [φ(x), Π(y)] = iδ(x−y) and φ has the meaning of a phase. The Hamiltonian is brought to the standard field theory notation [4] by a simultaneous rescaling of the fields φ → gφ and Π → g −1 Π, but in our context Eq. (1) is more convenient. Above, c is the light-cone velocity of correlation spreading, m tunes the overall mass scale, and g is the interaction. The parameter g governs the scaling dimension [3] (in an RG sense) of cos φ: for g 2 > 8π the cosine term is irrelevant and the theory flows to the free boson conformal point. For g 2 < 8π, interactions become relevant. As a consequence, the mass parameter m has an anomalous dimension scaling as [m] = 1 − g 2 /(8π). We postpone a more detailed discussion of the field theory to Section III. Here, we present its realization in tunnelcoupled spin chains (see Fig. 1) following the method of Ref. [21]. We thus consider two weakly tunnel-coupled XXZ spin chains where (3) with α =↑, ↓ and (S x iα , S y iα , S z iα ) the spin-1/2 operators. We choose the coupling between the two chains H ⊥ in the form of a weak tunneling Here, the coupling J ⊥ is assumed to be small (specified below). To see the emergence of sine-Gordon physics, one first neglects the tunneling and focuses on the effective low-energy description of the two chains, assuming to be close to the ground state. To this end, one proceeds within the Luttinger Liquid approach [3,45,46] by introducing a phase field φ α and its conjugate field Π α for each of the two chains. It should be stressed that the Luttinger Liquid approach is of much wider applicability than the spin chain and only requires a U (1) conserved charge (the z−magnetization in this case) and gapless excitations. In order to have linearlydispersing gapless excitations over the ground state of Eq. (3), we focus on the XY-phase with ∆ ∈ [−1, 1].
Each spin chain is thus described by the Luttinger Liquid Hamiltonian where the Luttinger parameter K and sound velocity v s fully characterize the many-body interactions. Within the bosonization language, the spin operators are represented as [47] where only the most relevant terms (in an RG sense) are retained. The constant prefactor α in front of the phase is non-universal and depends on the microscopic properties of the model. We will treat it as a constant, that will however renormalize the bare mass of the sine-Gordon model, so it must be carefully taken into account. As a second step, one now reintroduces the coupling between the two chains within the bosonization approach H ⊥ = |α| 2 J ⊥ dx cos(φ ↑ − φ ↓ ). This strategy is only valid in the weak tunneling regime. Finally, one rotates the fields as φ ± = φ ↑ ±φ ↓ and Π ± = (Π ↑ ±Π ↓ )/2. In this new basis, the symmetric and anti-symmetric degrees of freedom are explicitly decoupled. While the former are still described by a non-interacting Luttinger liquid, the sine-Gordon Hamiltonian naturally emerges in the antisymmetric sector The equivalence with H SG (1) is readily established with the identification c = v s , cg 2 = 2π/K and J ⊥ |α| 2 = c 2 m 2 /g 2 . The sign of J ⊥ is unimportant, since it can be changed by a global shift φ − → φ − + π. Various tunnel coupled one-dimensional discrete or continuous models would equally serve our purposes for realizing the sine-Gordon model. However, considering XXZ spin chains has the great advantage that the Luttinger parameters are analytically available. Indeed, the XXZ spin chain is a well known integrable model whose exact solution in the zero magnetization sector gives [48] We tested the assumption of the decoupling between the symmetric and antisymmetric sectors by measuring correlation functions and numerically extracting the Luttinger parameter K of the symmetric sector in Appendix A. These estimates indicate that K in the weak coupling regime J ⊥ /J < 0.3, is very well compatible with the single-chain exact result, Eq. (8), showing that hybridization effects between the symmetric and antisymmetric sectors induced by beyond-sine-Gordon corrections are small. In the remainder of this paper we will at most set J ⊥ /J = 0.2 to be in the sine-Gordon regime and determine K and v s from Eq. (8).
We would like to briefly comment on the asymmetric choice in the interchain coupling (4): the manifest absence of a spin-spin interaction in the z-direction is not accidental. Indeed, such an additional term would spoil the sine-Gordon effective description, by coupling the symmetric and antisymmetric sectors in a highly nontrivial way [3,49]. While unproblematic on the theoretical ground, a complete suppression of the S z ↑ S z ↓ coupling can be challenging in experimental spin chain realizations. To this end, a more convenient experimental setup can be realized with the Bose-Hubbard model with on-site interactions and nearest-neighbor hoppings [44]. However, for our numerical studies the restricted local Hilbert space of the spin chain and the analytic solutions for its Luttinger parameters are advantageous over directly simulating the Bose-Hubbard model, which is why we focus on the former.

III. THE QUANTUM SINE-GORDON MODEL
We briefly review the main features of the sine-Gordon field theory and its exact solution, for an exhaustive dis-cussion see Refs. [4,5]. We mostly refer to the notation of Eq. (1), but occasionally relate g and K explicitly when convenient. As already mentioned in the introduction, the sine-Gordon field theory belongs to the class of exactly solvable models: its integrability is established both on the classical [50] and the quantum level [2]. The key characteristic of integrable systems is the presence of infinitely many local conservation laws [51]: this has deep consequences on the excitation spectrum and dynamics. First, conservation laws ensure the existence of asymptotic multiparticle states, in spite of the strong interactions. In the sine-Gordon case, the fundamental excitations are solitons and antisolitons connecting between the valleys of the cosine potential. These excitations have equal masses M s and relativistic energy E s (θ) = c 2 M s cosh θ and momentum P s (θ) = cM s sinh θ, with θ the relativistic rapidity. The soliton mass scale has a highly non-trivial dependence on the interaction and bare mass scale m [52] where ξ = (8π/g 2 −1) −1 = (4K−1) −1 , and Γ is the Euler-Gamma function. In addition to solitons, the model also features non-topological excitations called breathers. Their masses are quantized according to where n = {1, ..., N } with N = ξ −1 = 4K − 1 . Notice that the interaction g, or equivalently the Luttinger parameter K, tunes the quantumness of the model: for large values of K, the mass gaps between breathers diminish and ultimately merge in a continuum, which corresponds to the excitations of the classical theory. Hence, quantum effects are most prominent for small values of K (i.e. strong interactions in the one-dimensional channels).
Interactions among sine-Gordon excitations are encoded in highly non-trivial scattering matrices. Importantly, the conservation of infinitely many charges allows only for elastic scattering. In this respect, two excitations with different masses (e.g. two breathers) can only be transmitted and never reflected: in this case, the logarithmic derivative of the scattering matrix is readily connected with the post-scattering space displacement experienced by two colliding wave packets. In Section V we will explicitly study such scattering events. The twoparticle scattering matrix is the backbone of the integrability approach to the model and it is exactly known [2]: for the sake of completeness, we report the lengthy expression and the space-displacement interpretation in Appendix B.
The important task of connecting the asymptotic states with actual observables requires the knowledge of the matrix elements, also known as form factors. The form factor boostrap [5] exploits the analytical properties of the scattering data to compute the sought-after matrix elements: this procedure is extremely challenging, but when possible leads to exact results. For our purposes, we are mostly interested in the vertex operator e iφ− [53][54][55]. Here we focus on the matrix elements between the ground state |0 and one-breather excitations |B n (θ) (12) and the ground state expectation value G g = 0|e iφ− |0 is exactly known [56] In the next section, we begin to explore the low-energy phase space of the spin ladder and quantitatively compare it with the sine-Gordon predictions.

IV. SPECTRAL ANALYSIS OF THE LOW ENERGY SECTOR
We systematically probe the low-energy sector of the spin ladder numerically and aim for a quantitative comparison with the sine-Gordon results. This allows us to provide bounds on the validity of the effective low-energy description. As the main probe of the excitation spectrum we consider the spectral function where O j is a local operator centered on site j, |0 is the ground state of the model (2). The Gaussian en- smoothens out the Gibbs phenomenon in the Fourier transformation caused by the finite-time cutoff present in numerical simulations. To subtract the zero-frequency signal of the ground state, we consider the connected correlation function . In order to probe the particle content of the theory, we evaluate the spectral function for the single-rung op- . The nonuniversal prefactor α is the same governing the bare mass scale (7) and does not depend on the chain coupling J ⊥ . The spectral function (14) can be computed by inserting a Lehmann decomposition of the identity in between the two observables and by expanding in the basis of the asymptotic states of the field theory. By retaining only the one-particle contributions, one obtains As a first comparison, we track the relativistic dispersion law of the breathers and the mass ratios. To this end, we access the microscopic spectral function with matrix-product-state (MPS) techniques to evaluate the real-time correlation function [57,58]: first, we obtain the ground state |0 by the density-matrixrenormalization-group (DMRG) algorithm [42,59], then we perform a real-time evolution of the state O 0 |0 with the time-dependent-block-decimation (TEBD) algorithm [57,60,61]. For further details about the numerical evaluation of the spectral function, we refer to the Appendix C.
In Fig. 2 (a) we show the spectral function of a ladder with ∆ = −0.4, corresponding to K = 1.35, and J ⊥ = 0.2J. For this choice of parameters, we expect N = 4 breathers to appear in the spectral function as can be indeed observed. In Fig. 3 further spectral functions for a different number of breathers N = 2, 3 are shown and the coupling strength J ⊥ /J = 0.1, 0.2 is varied. Notice that, at least for small momenta, the signal nicely follows the expected behavior of the relativistic dispersion law E Bn (k) = c 4 M 2 Bn + c 2 k 2 . The feature observed around k π/2 is a clear deviation from the sine-Gordon description and is a remnant of the single-chain spinon dispersion. In the single chain the lower edge of the spinon continuum at k = π/2 is perturbatively given by (1 + 2∆/π)J [62,63]. The weak tunnel coupling is only expected to strongly modify the low-energy regime ω J. In the higher energy regime our correlator perturbatively contains terms of the form 0 ↑ |S − j↑ (t)S + 0↑ |0 ↑ with |0 ↑ the ground state restricted to the single chain, and unveils the lower edge of the spinon continuum at high momenta.
To test the validity of the mass prediction Eq. (10), we then extract the breathers energies from the peaks of the spectral function at k = 0. In this case, we use the soliton energy and K as fitting parameters; see Fig. 2 (a) inset. As best fit parameters we obtain K = 1.35 ± 0.03 and c 2 M s /J = 0.56 ± 0.01. The best-fit K is thus compatible with the single-chain exact result of K = 1.35 given by Eq. (8). The excellent agreement of numerical data with the mass law is already a clear signature of the sine-Gordon effective theory, but we now want to quantitatively compare the whole spectral functions Eqs. (14) and (15), thus confirming that breathers contribute through the form factor expression (11).
The sine-Gordon model determines the spectral function up to an overall constant α coming from the singlechain bosonization (6) which, in turn, sets the sine-Gordon bare mass scale and eventually the soliton mass through Eq. (9). Reversing the logic, we can look at the physical soliton mass as the free parameter to be determined, which unambiguously fixes the microscopic non-universal proportionality factor α. We first check the non-trivial dependence of the mass gap (mass of the lightest breather) with the coupling across the rungs J ⊥ . Indeed, by feeding the microscopic value of the bare mass in Eq. (9), we obtain ∆E = c 2 M B1 ∝ J (1+ξ)/2 ⊥ , where the proportionality constants are J ⊥ independent. In Fig. 4, we show the numerically computed energy gaps as a function of the rung coupling for several interaction strengths ∆, finding an excellent agreement with the non-trivial power law behavior in the whole range of explored parameters.
Equipped with this benchmark of Eq. (9), we can now pursue the program outlined above and extract the parameter α from the measurement of the energy gap. In Fig. 2 (b) we quantitatively compare the energy dependence of the spectral function (14) with the sine-Gordon prediction (15) for different values of the momenta k: at small momenta, the agreement is excellent. We stress once again that α is the only fitting parameter. The Gaussian broadening of the peaks governed by σ in Eqs. (14) and (15) has been set equal in both cases. As one moves to higher momenta, small deviations from the field theory prediction become evident, but the peaks of the excitations are still clearly seen (see e.g. k = 0.3π). In this regime, we interpret the deviation as lattice corrections to the sine-Gordon relativistic dispersion relation, which is slightly bend towards smaller energies. As the momentum is further increased, the spinon contribution of a single chain starts to become dominant and the energy peaks of the sine-Gordon breathers disappear.
We further analyze the spectral function for different parameters, changing the number of breathers in the sine-Gordon spectrum, in Fig. 3. To quantitatively probe the field theory, we show the comparison between the numerical and analytic spectral function at zero momentum k = 0. Reducing the value of J ⊥ , the energy scale of the excitations is lowered and the spectral function develops a continuous tail associated with multiparticle contributions, to be contrasted with the isolated peaks of single-breather terms. To capture this continuum, we add the soliton-antisoliton contribution to the spectral function (15). The solitonic form factors of the vertex operators have been computed in Ref. [53]. In Fig. 3 we compare the spectral function at k = 0 obtained by considering only the contribution of breathers (dashed line) with the extension to the soliton-antisoliton corrections (solid line) and the numerics (circles). The comparison shows that solitons capture the tails of the spectral function well. Depending on the value of K, the energy of the soliton-antisoliton pair can be comparable to the excitations of two breathers B 1 : we checked this correction to the spectral function, but the B 1 − B 1 contribution turned out to be negligible. We thus attribute the small discrepancy with the numerical data to corrections beyond sine-Gordon.
From the spectral functions (shown in Figs. 2 and 3) it becomes apparent that the spectral weight of the higher lying breathers, with breather index n > 2 is strongly suppressed with respect to the lower lying ones. In Appendix D we provide a form factor argument to increase the relative weight of the higher lying breather, that is based on computing the spectral function over a correlation function of a linear combination of single-and tworung operators.
While the spectral function (14) gives direct access to the low energy spectrum of the coupled chains, measuring unequal-time correlation functions is challenging in realistic experimental scenarios. To this end, the fluctuationdissipation theorem can be used to relate the spectral function to the energy absorption due to periodic drives. We consider a periodic modulation of the Hamiltonian H(t) = H + V (t), where H is given by Eq. (2), and where we chose the perturbation as a linear combination of two terms . The choice of the drive is motivated on physical grounds. O 1 is the z−magnetization difference between the two chains. Hence, driving O 1 can be achieved by applying a small magnetic field of opposite sign to the two chains. On the other hand, modulations coupled to O 2 are obtained by manipulating the interchain coupling (4). Here, we focus on the linear response on the top of the ground state. Hence we limit ourselves to the case < J ⊥ . Within bosonization (6), O 1 ∼ dx Π − (x) and O 2 ∼ dx cos(φ − ). Therefore, the operators O 1 and O 2 respectively excite odd and even breathers from the ground state. In linear response theory, the fluctuationdissipation theorem links the averaged absorbed energy to the spectral function at zero momentum. Defining where T = 2π/ω is the driving period. Hence the spectral function is directly probed by the slope of the on-average linear growth of the energy. Notice that the different parity of the operators O 1,2 leads to a splitting of the spectral function S(ω, k) = a 2 S O1 (ω, k) where S Oµ is the spectral function of the operator O µ . Hence, the choice of the coefficients a, b does not have any physical meaning within the linear regime and we consider a = b = 1 to probe the even and odd sectors within a single measurement. In Fig. 5 (a) we show the energy growth for various driving frequencies, and in Fig. 5 (b) we show the corresponding response function extracted from linear fits in several time windows. We compare the results to the breather energies extracted from the two-point correlation function, illustrating a very good agreement. The spectral weight of B 3 and B 4 is probably too low to render these excitations observable in experiments, however the ratio between B 2 and B 1 is already sufficient to obtain the mass scaling.

V. REAL-TIME SCATTERING OF SINE-GORDON EXCITATIONS
The excellent comparison between the ladder spectral function and the sine-Gordon field theory allowed us to identify the regime of validity of the latter. As a next step, we probe genuine nonequilibrium dynamics of the system. A clear signature of integrability can be found in scattering events: the existence of infinitely many conservation laws ensures complete elastic scatterings even if non-elastic scatterings are, in principle, energetically allowed. We thus create localized wave packets of the lighter breathers and observe their scattering dynamics. In this case, the sine-Gordon scattering matrix ensures complete transimission of breather-breather scatterings (to be contrasted with soliton-antisoliton scattering which, due to additional symmetries, can also be reflective). Complete transmissive scattering arises also in non-interacting theories, but in the sine-Gordon field theory interactions are strong and manifest as a spatial displacement after scattering. The displacement is completely determined by the exactly known scattering matrix [2], as we briefly report in Appendix B (see also Refs. [64,65] for numerical computation of scattering displacement in integrable spin chains).
The sine-Gordon excitations are collective low-energy modes that are not trivially connected to the spin configurations. To engineer a local breather excitation, we create a Gaussian wave packet with momentum k, centered around site d with widthσ by applying the operator on the ground state O (1) d,k |0 . Within bosonization S z j↑ − S z j↓ Π − , hence the low-energy excitations are determined by the sine-Gordon form factors B n (θ)|Π − |0 which are exactly known [5], but their analytical expression is not needed for our purposes. The Π − operator is odd under the Z 2 symmetry of the sine-Gordon model φ − → −φ − . Moreover, the matrix element is non vanishing only for breathers odd under parity, which correspond to odd integers n, and gets smaller as n is increased. Therefore, the operator (18) will with highest probability couple to the lightest breather. But this is not the only contribution: heavier breathers are expected to be excited as well and, most importantly, high-energy excitations beyond the sine-Gordon description are inevitably produced. To eliminate the unwanted contributions and target low-energy excitations within an energy interval ∆E, we perform a short imaginary-time evolution e −τ H , τ ∼ 1/∆E. After this, we indeed obtain clean wave packets that are moving with a speed compatible with the group velocity v 1 (k) = ∂ k E B1 (k). For further details about the creation of the wave packets, we refer to Appendix E.
With the operator O d,k Eq. (18), we can target the first breather B 1 , but we also want to scatter wave packets of different species to highlight the sought-after transmission. To this end, we also create wave packets by using another operator O In this case, we need to explicitly impose orthogonality with respect to the ground state 0|O (2) d,k |0 = 0 by subtracting the expectation value o j = 0|S + j↑ S − j↓ + h.c.|0 . In bosonization, [S + j↑ S − j↓ +h.c.] ∼ cos φ − , which is an even operator. Hence, in the low-energy sector O (2) d,k |0 couples only to wave packets of even breathers. Therefore, by cooling down the state O scattering wave packets in the upper-half (x u-h ) and in the lower-half (x l-h ) of the system as obtained by Note that these measures are only valid in the asymptotic regions before and after the scattering event, where the separation between the wave packets is much larger than their typical interaction length. Comparing the average position between the interacting and freely propagating case, we find a displacement of δx 0.5 (see inset in Fig. 6(c)). This value can be compared to the analytical result obtained from the S-matrix of the sine-Gordon theory. To obtain a concrete result, we have used the following: (i) the Luttinger parameter K and the sound velocity v s = c as obtained from the single chain by Eq. (8), for ∆ = −0.4, K = 1.35 and c = 0.73J, (ii) the soliton mass M s = ∆E/(2c 2 sin(πξ/2)) where we determined the energy gap ∆E/J = 0.26 from the spectral function, (iii) the average speed of the wave packets v = ±0.60J which we have obtained from a fit of the freely propagating wave packets and not from the original input momentum chosen in Eq. (18) as the imaginary time evolution could potentially slightly lower the initial momentum of the wave packet (see Appendix E). Finally, when the two wave packets collide, each excitation will scatter with all the breathers contained in the other wave packet and pile up the cumulative space displacement. Hence, we need to estimate the number of excited breathers. To this end, we consider the total energy carried by each wave packet and divide by the breather dispersion law computed at the rapidity extracted from the velocity measurement. From this analysis, we found that we excite approximately one breather in each wave packet. Using this last piece of information, we then obtain an analytic value of δx d,k1 |0 with k 1 = −0.2π, k 2 = 0.15π and d = 12. The wave packet associated with B 2 (the heavier quasi-particle) is expected to have a smaller velocity than the B 1 wave packet, and therefore also a larger broadening in time. The energy as a function of time for this process is shown in Fig. 6(d). In Fig. 6(f), we plot the average position of the wave packets. As the wave packet associated with B 2 is moving slower we extend (decrease) the upper-half (lowerhalf) of the system by 10 sites in the measurement of Eq. (20) (u-h/l-h = 60/40), but only after the scattering. Before the scattering we still divide the system in halves (u-h/l-h = 50/50). To check that is indeed a valid choice, we have also fitted a bi-modal Gaussian distribution to our data before and after scattering (no such fit is possible during the scattering event). The mean values of this distribution are in good agreement with the positions obtained from Eq. (20). From our data we extract the displacements δx B2 = 1.18 again in very good agreement with our numerical estimates. In this case, we have also checked explicitly that the asymptotic velocities of the wave packets after scattering are not compatible with a reflection event, i.e. with the non-transmissive solution of the energy-momentum conservation E 1 (k 1 )+E 2 (k 2 ) = E 1 (k 1 ) + E 2 (k 2 ) and k 1 + k 2 = k 1 + k 2 . Indeed, plugging in the above-mentioned velocities would give v 1 = 0.55J and v 2 = −0.45J as velocities after reflection. This is in contrast to the B 1 − B 1 event, for which reflection and transmission are degenerate.

VI. CONCLUSIONS & OUTLOOK
The sine-Gordon field theory is a ubiquitous lowenergy approximation of one-dimensional quantum systems. Recent experiments focused on implementing it on highly-controllable quantum simulators based on coupled one-dimensional quasi-condensates [21,23]. In our work, we considered two tunnel-coupled spin chains, realizing a lattice version of the sine-Gordon model proposed in Ref. [21]. The lattice realization has advantages compared with the quasi-condensates: first, the one-dimensional chains can be strongly interacting, leading to a sine-Gordon model deep in the quantum regime in contrast to the weakly interacting quasi-condensates. Second, the dynamics of coupled chains can be accurately determined numerically, which allows for a direct observation of the emergent sine-Gordon description and for a quantitative analysis of its validity regime.
By a numerical characterization of the spectral function, we have studied the sine-Gordon spectrum in the low-energy sector of the chain. Comparing the numerical data with exact predictions, we obtain precise regimes for the validity of the field theory. The observation of the rich sine-Gordon spectrum already supports the field theoretic approximation. Yet, a striking feature of the model is its integrability, resulting in unconventional nonequilibrium properties. As a probe of integrable dynamics, we studied scattering events of low-energy wave packets. Not only we observe fully transmissive scattering as expected from integrability, but we also indirectly observed the sine-Gordon scattering matrix by measuring the wave packet displacement after scattering.
Our work demonstrates that coupled ladder geometries are useful to realize a simulator of the sine-Gordon model in the quantum regime. Similar realizations as the one considered here based on coupled Bose-Hubbard models can be realized with current experimental techniques [44], leading to various interesting questions for future investigations. First, a characterization of experimentally-feasible protocols to selectively create and manipulate the sine-Gordon excitations will be of interest. To this end, numerical benchmarks will play an important role to assessing the validity of the quantum simulation and to keep unwanted effects under control. Second, the important question of characterizing the sine-Gordon dynamics beyond the few-excitation regime presents a formidable challenge. To this end, the theory of Generalized Hydrodynamics developed for integrable models [66][67][68] promises access to large scale dynamics. Numerical benchmarks will be helpful in testing the hydrodynamics predictions on small scales and in characterizing the short-time transients that are beyond the hydrodynamic approximation.

VII. ACKNOWLEDGEMENT
We are grateful to F. Essler for helpful remarks on the bosonization of the coupled spin chains, and thank I. Lovas, L. Vanderstraeten and W. Kadow for useful discussions. Tensor-network calculations were performed using the TeNPy Library [43]. We acknowledge support from the Deutsche Forschungsgemeinschaft As discussed in Section II, the bosonization of the coupled chains predicts the emergence of the sine-Gordon field theory in the odd degrees of freedom, while the even sector remains a gapless Luttinger liquid. In this section, we provide a sanity check of this claim by numerically probing the static correlation functions of the ground state and comparing them with their analytic forms given in Refs. [69,70]. Indeed, since the contribution of the gapless sector dominates the exponentiallydecaying contribution of the massive antisymmetric part φ − , the critical behavior is already clearly seen in singlechain correlation functions, furthermore giving access to K. The results are shown in Fig. 7. These fits are however hard to perform, particularly in the regime where ∆ > 0 due to the large oscillations in the correlation functions, which results in the large errorbars. We therefore show the results for both and with r = |j − k|. Notice that we include the interchain correlations as well, although we are considering the weak coupling limit. We have found that by including them the fits were improved, particularly for zz r . The errorbars shown in the figure are a combination of the uncertainty on the fit, and of the averaging over some different fitting domains [r min , r max ]. We incorporate the latter because the fit of zz r is quite sensitive to the choice of r min , due to non-universal short ranged corrections.
These estimates confirm that the Luttinger parameter of our model agrees with the single chain exact result (8) in the weak tunnelling regime. In particular, deviations are only found around J ⊥ /J 0.3. For our analysis in the main text, we have considered J ⊥ /J = 0.1, 0.2.

Appendix B: The sine-Gordon scattering matrix and wave-packet displacement
We report the sine-Gordon scattering matrix and its connection with the spatial displacement after a scattering event has occurred. The exact expression of the scattering matrix can be found in Ref. [2]. From that the breather-breather scattering matrix can be obtained by the form factor bootstrap. Since we are only interested in breather-breather scattering events, we focus on the latter for simplicity. Let S (n,m) (θ) be the scattering matrix of a scattering between a breather of species n and where the renormalized interaction ξ has been defined below Eq. (9). The case n < m can be recovered by using the unitary relation S (n,m) (θ)S (m,n) (−θ) = 1. The scattering matrix is a complex number with modulus 1, therefore one can define the scattering phase Θ (n,m) (θ) = −i log S (n,m) (θ). The derivative of the scattering phase is connected with the spatial displacement after scattering, see e.g. Ref. [71]. In particular, let δx (n,m) Bn (θ n , θ m ) be the spatial displacement experienced by a breather of species n with rapidity θ n colliding with a breather of species m and rapidity θ m , then δx (n,m) Bn (B2) This is the analytical expression of the spatial displacement used in Section V. the Fourier transformations to momentum and frequency space. To obtain the time evolution, we apply the TEBD algorithm [57,60,61] to the state O 0 |0 , were |0 is the ground state of the ladder (2) obtained from DMRG [42,59]. We group the sites connecting each rung of the ladder (labelled by ↑, ↓) to a single site, and map the ladder system to a chain with an increased local Hilbert space of dimension d = 4. Within this chain there are only two-site interactions. Hence, we can apply a fourth order TEBD scheme efficiently. As a next step, it suffices to compute the overlaps with O j |0 at different time steps as with E 0 the ground-state energy. The operator O 0 acts in the middle of the system and ideally we want to exclude boundary reflections. Therefore, the maximal time of the simulation is chosen as t max L/(2v s ) with v s the sound velocity. As the resolution in frequency space is directly related to t max , we simulate the largest system sizes possible. In practice, we simulated ladders up to L = 151 rungs for moderate ∆ ∈ [−0.4, 0.1]. The sound velocity is set by ∆ (see Eq. (8) in the main text), so also t max and the frequency resolution will depend on ∆.
As the sound velocity becomes larger with increasing ∆, the entanglement growth will be also more pronounced.
In Fig. 8 we show the entanglement growth in a chain of L = 151 grouped sites, with J ⊥ /J = 0.1, for two different values of ∆ each at three different bond dimensions, where (C2) with b running over all bonds, and ρ b the bipartite reduced density matrix with respect to that bond. In the simulations shown in the main text, we typically use a bond dimensions of D = 1500. Before taking the Fourier transform to frequency space, we multiply the data with a Gaussian window with standard deviation σ [57], which we tune as function of the maximal time of the simulation t max . This sets the resolution in frequency space, and induces a broadening of the peaks.

Appendix D: Form factor suppression of B1
The relative spectral weight of the breathers in the spectral function can be tuned by considering a linear combination of a single and a double rung operator, instead of just the single rung operator. This could be useful because the spectral weight of the breathers B n with breather index n > 2 is strongly suppressed compared to the lower lying ones B 1 , B 2 . Here we discuss a sine-Gordon form factor argument to reduce the weight of B 1 compared to B 3 (we will work in the odd breather sector). Introducing the single and two rung operators which according to bosonization are (α 1,2 ∈ R) Our goal is to find a linear combinationÕ = O 1 + AO 2 of these operators such that Within sine-Gordon, we have the following form factor for B 1 − 0 We have that and therefore to satisfy Eq. (D3) A − 1 2 cos(πξ) Due to the interchain symmetry 0|O 1,2 |0 = 0|O 1,2 |0 * , we can simply translate the argument to the anti-symmetric combination of the operators Such that when we compute the spectral function Eq. (14) with O j (t) =Õ j (t), B 2 vanishes because of parity, and B 1 is significantly reduced by the above construction. We have numerically computed the four contributions to the correlation function separately Their respective spectral functions at k = 0 are shown in Fig. 9 for different values of ∆. Each of the spectral functions detects B 1 and B 3 , however, with different relative weights. In the spectral function obtained from the two-rung operators, Fig. 9(d), there is also a continuum contribution present. When summing up the correlation functions according to Eq. (D9) we can indeed reduce the contribution of the first peak significantly. The resulting spectral function at k = 0 are shown in Fig. 9  The action of the operators given in Eqs. (18) and (19) on the ground state does not only excite the low-energy sine-Gordon part of the spectrum, but also higher-energy contributions. Indeed, in Fig. 10(a), we show the time 1J. In the case without imaginary time evolution, there is a two-particle contribution beyond sine-Gordon, that moves with a velocity close to the sound velocity vs (pink line). A short imaginary time evolution eliminates this high-energy contribution, and we recover a single mode which moves with a velocity compatible with the group velocity as read off from the B2 dispersion (orange line).

evolution of the state O
(2) 0,k2 |0 with k 2 = 0.15π and |0 the ground state. There is a highly-energetic two-particle contribution, as there is also an outgoing signal with negative momentum, that blurs the signal of the sine-Gordon breather. To eliminate this high-energy continuum contribution, we 'cool down' the system by performing a short imaginary-time evolution e −τ H , maximally up to imaginary times of around τ ≈ 4/J. After this, we indeed obtain a clear single mode that is moving with a speed lower than the sound velocity and compatible with the group velocity v(k) = ∂ k E B2 (k) according to the massive-relativistic dispersion of the second breather, as illustrated in Fig. 10(b).
We estimate the group velocities of our wave packets from fitting the massive-relativistic dispersion relation to the individual breather dispersions as obtained from the spectral functions. We restrict those fits to small momenta. We find that the fitted c of B 1 and B 2 slightly differs, as shown in Fig. 11(a). To enforce consistency, we fixed c to the single chain result of the sound velocity v s , given by Eq. (8). For ∆ = −0.4, this corresponds to v s = 0.726J. This value is also approximately the average of the fitted c's shown in the legend of Fig. 11(a). We obtain the breather mass from the fitted value of the gap via M Bn = ω Bn (0)/c 2 , as indicated by the green lines.
There are no conserved quantities during the imaginary time evolution to cool down the state (this is always the case with TEBD-based algorithms). In addition, as imaginary time evolution is non-unitary the energy will not be conserved. In the long-time limit, imaginary time evolution acts as a projector on the ground state. How- ever, we are not interested in the long time limit, and start from a state that is orthogonal to the ground state. By choosing the imaginary time-step small enough, we used dτ = 0.0025J in our simulations, the state stays orthogonal to the ground state (at least up to the relevant time scales). After the short imaginary time evolution, we checked that (i) the obtained state is still orthogonal to the ground state, and (ii) that the energy E(τ ) of the obtained state is roughly compatible with E(τ ) − E 0 ≈ ω B1 (k 1 ) + ω B2 (k 2 ) with E 0 the ground state energy. The total energy difference as a function of imaginary time is shown in Fig. 11(b). The imaginary time evolution can in principle also lower the energy of the excitation along the single-particle breather dispersion, such that the group velocity of the wave packet could be overestimated when obtained from the initial momenta that are put in Eqs. (18), (19). We have checked this for our data, and this effect seems to be present but not very strong. When a higher precision is required, we directly fitted the velocity of our wave packets instead of determining it via the input momentum.