The quantum sine-Gordon model with quantum circuits

Analog quantum simulation has the potential to be an indispensable technique in the investigation of complex quantum systems. In this work, we numerically investigate a one-dimensional, faithful, analog, quantum electronic circuit simulator built out of Josephson junctions for one of the paradigmatic models of an integrable quantum field theory: the quantum sine-Gordon (qSG) model in 1+1 space-time dimensions. We analyze the lattice model using the density matrix renormalization group technique and benchmark our numerical results with existing Bethe ansatz computations. Furthermore, we perform analytical form-factor calculations for the two-point correlation function of vertex operators, which closely agree with our numerical computations. Finally, we compute the entanglement spectrum of the qSG model. We compare our results with those obtained using the integrable lattice-regularization based on the quantum XYZ chain and show that the quantum circuit model is less susceptible to corrections to scaling compared to the XYZ chain. We provide numerical evidence that the parameters required to realize the qSG model are accessible with modern-day superconducting circuit technology, thus providing additional credence towards the viability of the latter platform for simulating strongly interacting quantum field theories.


Introduction
The investigation of strongly interacting complex quantum systems remains one of the outstanding challenges of modern physics. Despite the remarkable progress on both numerical as well as analytical fronts, systematic and well-controlled non-perturbative analysis of many quantities of interest remain intractable. Quantum simulation provides a promising alternative to the aforementioned conventional techniques towards tackling these problems [1]. There are two approaches to quantum simulation: digital and analog. In principle, digital quantum simulation is universal [2,3]. It can be performed by a digital quantum computer built out of qubits. A digital quantum simulation of a many-body Hamiltonian comprises encoding of the target Hamiltonian as a trotterized sequence of one and two-qubit gates and readout of desired observables. However, opening up the digital quantum computer to externally applied gates and measurement apparatus leads to unavoidable decoherence of the physical qubits constituting the quantum computer. To combat for the finite lifetimes of the physical qubits as well as imperfections of the applied gateset, quantum error-correction [4,5] is essential. In the recent years, spectacular process has been made towards realizing such a universal computing machine. In fact, noisy, intermediate-scale, quantum machines have already been shown to be capable of simulating certain aspects of fewbody systems [6,7]. However, extrapolating these efforts to the many-body domain will require enormous overhead and is likely to remain elusive in the immediate future. A more near-term approach to simulate many-body physics is analog quantum simulation, where a given quantum system is tailored to simulate another [8,9,10,11,12,13]. In this approach the target Hamiltonian is realized by specifically engineering the given quantum system and letting it naturally evolve with time. There is no need to implement trotterized gate-sets. Thus, the physical degrees of freedom that comprise analog simulators do not need to be individually addressed and can be isolated from losses much more than their counterparts in a digital quantum computer. The imperfections in the engineered quantum system lead to 'errors' which are typical of an experimental realization, e.g., finite correlation lengths at criticality due to the disorder, etc. One of the biggest advantages of analog quantum simulation is the availability of a wide range of viable experimental platforms. Analog simulators based on trapped atoms have been used to experimentally simulate strongly correlated systems, topological phases of matter and gauge theories [14,15,16,17], while trapped-ion based simulators have investigated problems in quantum magnetism [18,19]. Finally, superconducting quantum electronic circuit (QEC) based simulators have been used to experimentally probe quantum electrodynamics in the strong and ultra-strong coupling regimes [20,21,22,23].
In this work, we advance the research direction of faithful analog simulation of quantum field theories (QFTs) with QEC lattices [13]. Here, faithful refers to the fact that the degrees of freedom of the QFT are faithfully represented by the underlying lattice degrees of freedom and do not arise out of mathematical manipulations like bosonization. This is particularly important in multi-field QFTs, where properties as fundamental to a QFT as integrability, can be difficult to relate in the fermionic and the bosonic counterparts [24,25,26,27,28]. With this potential generalization to multi-field situations in mind, here we investigate, with the density matrix renormalization group (DMRG) technique, a faithful analog QEC simulator for one of the paradigmatic integrable QFT models: the quantum sine-Gordon (qSG) model in 1+1 space-time dimensions. The QEC simulator is a one-dimensional array of suitably-arranged Josephson junctions and provides a faithful lattice-regularization for the qSG model using only local interactions. In the first part of the paper, we provide numerical evidence that the long-wavelength properties of the QEC lattice model are indeed described by the qSG field theory by computing various zero-temperature thermodynamic properties of the lattice model and comparing with analytical field-theory predictions. In the second part of the paper, we analyze the entanglement spectrum of the qSG model. To that end, we first compute the entanglement spectrum using the integrable Luther-Lukyanov lattice regularization [29,30] involving the quantum XYZ spin-chain [31]. The entanglement spectrum of the XYZ chain is related to the spectrum of the corner-transfer-matrices (CTMs) of Baxter's eight-vertex model [31] and thus, can be computed analytically [32,33]. Then, we calculate the entanglement spectrum of the QEC incarnation of the qSG model using DMRG. We show that the entanglement spectra of both the XYZ and the QEC regularizations comprise equidistant levels. We argue that the QEC and the XYZ entanglement level spacings are linearly related and verify this claim with numerical predictions.
The analyzed QEC simulator provides a different lattice-regularization of the qSG model compared to the known existing ones. The primary motivation behind analyzing this lattice model is its experimental feasibility. It is a straightforward generalization of the current experimental works which have so far realized the free, compactified boson conformal field theory (CFT) [20,21,22,23]. In fact, QEC systems which realize the qSG model in the semi-classcial limit have also been fabricated and experimentally analyzed [34,35,36]. However, as is shown in this work, the QEC regularization of the qSG model is also of intrinsic theoretical interest. In contrast to the XYZ chain regularization, where the qSG model arises out of Jordan-Wigner and bosonization transformations, in the QEC incarnation, the compact bosonic field φ, is faithfully represented at the lattice level. This makes the QEC regularization more easily generalizable to multi-field scenarios [13]. At the same time, the correlation functions of vertex operators e iβφ , as computed using the QEC regularization, are less susceptible to corrections to scaling, where β is the qSG coupling constant. These corrections arise when the correlation length of the underlying lattice model is not large enough compared to the lattice spacing; later we provide a quantitative estimate of when this effect becomes noticeable. The QEC model is more immune to the corrections to scaling compared to the XYZ chain since we start directly with bosonic fields on the lattice. This is in contrast to the XYZ chain, where spin-operator σ + is proportional to the qSG vertex operator e iβφ/2 to leading order [30,37,38]. The situation for the XYZ model is further worsened by the fact that the qSG coupling and the mass of the soliton cannot be independently tuned unlike in the QEC lattice model. Additionally, we show that the same issues concerning reaching the scaling limit in the XYZ model also plague the entanglement spectrum -for certain parameters, the XYZ model does not provide meaningful predictions for the qSG model, in contrast to the QEC model. Note that there is one other well-known regularization, proposed by Bogoliubov, Izergin and Korepin, of the qSG model based on nonlocal interactions of the bosonic degrees of freedom [39,40]. While the latter model is precious for analytical computations using the quantum inverse scattering method, its direct implementation in physical systems is not straightforward. Lastly, lattice regularizations based on QECs are amenable to DMRG and thus, provide a crucial tool to investigate properties of QFTs (e.g., entanglement between spatially separated regions, etc) which are not easily accessible using alternative methods like the truncated conformal space approach [41,42,43].
The article is organized as follows. In Sec. 2, we briefly summarize the relevant properties of the qSG model. In Sec. 3.1, we describe the QEC lattice model and provide approximate expressions for the qSG parameters in terms of the lattice parameters. In Sec. 3.2, we briefly recount the lattice-regularization based on the XYZ spin-chain. In Sec. 4, we provide DMRG and analytical results for the one-point and two-point correlation functions. Secs. 4.1 and 4.2 describe the DMRG results for the QEC and XYZ regularizations respectively. The effects of corrections to scaling for the two models are discussed in Sec. 4.3. Finally, in Sec. 5, we compute the entanglement properties of the qSG model using DMRG. We first present the analytical and DMRG results for the XYZ model in Sec. 5.1, followed by the DMRG results for the QEC lattice in Sec. 5.2. The consequences of the corrections to scaling on the entanglement spectrum for the two models are presented in Sec. 5.3. Sec. 6 presents a concluding summary and outlook. In Appendix A, we provide the analytical results of the zero-temperature two-point function for the qSG model.

The quantum sine-Gordon model
In this section, we briefly summarize the basic properties of the qSG field theory that will be relevant for this work. More details on various properties can be obtained in, e.g., Refs. [44,45,46,47,48,49].
The qSG field theory is an integrable deformation of the free, compactified boson CFT. Its euclidean action is given by where M 0 is the mass-parameter of the action and β is the coupling constant. We set = 1 throughout this work. We restrict ourselves to the regime when β 2 ∈ (0, 8π). In the classical limit, which corresponds to β → 0, the coupling constant β plays no role and can simply be scaled out. The resulting theory is the well-known classical sine-Gordon theory, which supports traveling wave-packet solutions, which propagate undistorted through the nonlinear wave-medium and scatter with only phase-shifts [50]. Our interest, however, is in the quantum regime, when the parameter β determines the spectrum of the theory. For β < √ 4π, the qSG model is in the attractive regime, where the spectrum of single-particle excitations consists of solitons, anti-solitons and breathers. For β > √ 4π, the model is in the repulsive regime when the spectrum of singleparticle excitations consists only of solitons and antisolitons. The fermionized version of this model corresponds to the massive Thirring model [44,45,51], where the solitons are the fermions in the Thirring model. The choice β = √ 4π corresponds to the free, massive (complex) fermion QFT.
Many exact results are available for the qSG model, including its exact spectrum and the Smatrix. The following ones are relevant for our work. The soliton mass can be derived to be [52] The mass of the n th breather state is given by The ground state energy density with respect to the free, compactified boson CFT is [52] E 0 (M) = − M 2 4 tan πξ SG 2 (4) and the ground-state expectation value of the local vertex operator e iβφ is given by [53] Note that the last two predictions are for the continuum theory and in general, hold for the corresponding operators on the lattice only up to proportionality constants that depend on β; see, for example, similar computations for the XXZ model [30] and the 2D classical Ising model in a magnetic field [54]. Furthermore, we have used in the above expressions, the standard CFT normalization (note the difference in definition of β compared to Ref. [53]):

Lattice regularizations for the quantum sine-Gordon model
In this section, we describe two lattice-regularizations for the qSG model. We start with the QEC lattice and discuss this in detail. Then, we briefly summarize the XYZ regularization following Ref. [30].

QEC model
The QEC lattice model for the qSG model comprises a 1D array of mesoscopic superconducting islands [see Fig. 1 (b)]. Two neighboring islands are separated by Josephson junctions (indicated by the green crosses within brown squares) with junction energy E J and charging energy E C J = 2e 2 /C J . In addition, each island is also separated from the ground plane by a Josephson junction (indicated by a purple cross within a red box), with junction energy E J 0 and charging energy E C 0 = 2e 2 /C 0 . For reference, we show the QEC lattice for the free, compactified boson CFT in Fig. 1 (a). For the latter, the Josephson junction on the vertical link in each unit cell is replaced by a capacitance C 0 . Throughout this work, we consider a homogeneous array with zero disorder in the off-set charges on the different superconducting islands. 1 The Hamiltonian describing the array is given by Here, the first term arises due to the finite charging energy of the mesoscopic islands and n i is the excess number of Cooper pairs on the i th island. 2 The finite junction capacitance C J leads The Josephson junction on the horizontal link leads to a β 2 ∈ (0, 8π) [see Eq. (13)] while the Josephson junction on the vertical link gives rise to the cosine potential of the sine-Gordon model. In both cases, the bosonic field is the continuum version of the discretized superconducting phase ϕ i (t) at the node i of the lattice. Physically, the i th node denotes the i th superconducting island.
to, in principle, infinite-range interaction between any two islands with a magnitude that decays exponentially with distance [55]. The relevant length scale is given by a √ C J /C 0 , where a is the lattice spacing. However, for realistic system parameters [56], it suffices to include only the nearest neighbor interaction [57], indicated by the second term in Eq. (7) with δ being a small parameter < 1. The third term arises due to the presence of a gate voltage E g at each node. The fourth term in Eq. (7) describes the coherent tunneling of Cooper-pairs between neighboring islands. The last term describes the tunneling of Cooper pairs across the Josephson junctions on the vertical links and is responsible for the qSG nonlinearity. Note that the operators n i , ϕ j are canonically conjugate satisfying [n i , e ±iϕ j ] = ±e ±iϕ j δ i j .
First, consider the case E J 0 = 0 [this corresponds to Fig. 1 (a)]. Then, the QEC lattice realizes a version of the Bose-Hubbard model with nearest-neighbor repulsion [57,58,59], where the role of the bosons is played by Cooper-pairs. Since the number of excess Cooper-pairs can be both positive and negative, the model reduces to a Bose-Hubbard model for quantum rotors, with the number of Cooper-pairs being conserved. The phase-diagram of this model has been analyzed perturbatively [57] and with DMRG [60]. For E J E C 0 , the system is in a Luttinger liquid (LL) phase of Cooper-pairs, where its long wavelength properties are well-described by the action: 3 Here, u is the plasmon velocity and K the Luttinger parameter. Since the lattice model is nonintegrable, exact expressions for u, K are not known in terms of the lattice parameters. Perturbative estimates for E J E C 0 are given by [55] u a 2E C 0 E J , K 1 2π Lowering E J /E C 0 causes the system to transition into either a Mott-insulating or a charge-densitywave phase, with pinned densities ρ = m/n. In the Mott-insulating phases, n = 1 and m ∈ Z, while for the charge-density-wave phases, n = 2 and m = 2k + 1, where k ∈ Z. The transition out of the LL phase to either of the other two phases with fixed density ρ 0 = m/n is caused by a perturbation of the form where θ is the field dual to ϕ and ρ is the particle (Cooper-pair) density on each superconducting island. From dimensional analysis, it follows that the transition at fixed density ρ = ρ 0 occurs at a Luttinger parameter of K c = n 2 /2, while that with a change of density occurs at K c = n 2 [61,62,63]. For the model with only nearest-neighbor interactions, n can be at most 2. Thus, in the LL phase, the Luttinger parameter K ≤ 4. It might appear surprising that we use Josephson junctions on the horizontal link, but consider only the limit E J E C J , E C 0 , when the nonlinearity of the Josephson potential in the fourth term of Eq. (7) plays no role. Indeed, we do not use the nonlinearity of the Josephson junctions on the horizontal links. However, the kinetic inductance associated with the Cooper-pairs, leads to a Luttinger parameter up to 4. If linear electromagnetic coil inductances were used instead of the Josephson junctions on the links of the array, apart from making the lattice boson non-compact, the Luttinger parameter would be ∼ Z/R Q ∼ 0, where Z is the impendance of the array (∼ 50Ω) and R Q is the impendance quantum (∼ kΩ) [56,64]. As explained below, having K ∼ 1 is crucial for exploring the quantum regime of the qSG model. Now, consider the case E J 0 0, which breaks the particle number conservation of the quantum rotor Bose-Hubbard model described above and contributes an additional term to the action: Rescaling the field ϕ and the space-time axis: we arrive at the action of the qSG model: where ϕ = βφ. Furthermore, the qSG coupling and the mass-parameter of the action are Note that the free-fermion point of the qSG model occurs at β = √ 4π, which corresponds to K = 4 and not K = 1. From the above considerations, thus, the QEC lattice model of Eq. (7) is limited to the attractive regime of the qSG model: β ≤ √ 4π since it has only nearest-neighbor interactions. However, in an actual experimental realization, the interaction can, in principle, be long-ranged, but exponentially decaying. The range of the interaction is determined by the ratio C 0 /C J [55]. The longer-range model, which supports larger Luttinger parameters, is important for realizing the repulsive regime of the qSG model and can be similarly analyzed generalizing this current work.
Note two important features of the QEC lattice model for the qSG model. First, the underlying lattice degrees of freedom accessible to numerical simulations are the vertex operators e iϕ i , which, to leading order, up to a proportionality constant depending on β, are the same as the QFT vertex operators e iβφ . Second, the deviations from the leading qSG action are captured by the term given in Eq. (10). They are irrelevant in the renormalization group sense, but more importantly, do not renormalize the relevant mass term in the qSG action. This allows us to have an analytical control on the mass-parameter of the continuum qSG action in terms of the QEC lattice parameters.

XYZ model
Now, we briefly summarize the qSG limit of the XYZ chain following Ref. [30]. Consider the Hamiltonian: where J x,y,z are the coupling constants and σ x,y,z are Pauli-matrices. We consider the parameter regime: J x ≥ J y ≥ |J z |. For J x = J y , the excitation spectrum is gapless and the low-energy properties of the model are described by a Luttinger liquid action [63]. The correlation length ξ XYZ diverges as where we denote the Luttinger parameter of the spin-chain by K XYZ . Close to the critical point, the system is described by the qSG field theory, with the action of Eq. (12). In this case, the qSG coupling and the mass of the soliton are given by where a is the lattice-spacing. For the continuum theory to be applicable, it is necessary to have ξ XYZ 1, which in turn restricts the achievable values of M. Note that for the XYZ chain, the free-fermion point of the qSG model occurs at K XYZ = 1. To leading order, the spin-creation operator σ ± (x) can be identified with the vertex-operator of the qSG model [30]: where the dots indicate corrections to scaling arising from irrelevant terms. More details on the corrections to scaling and their effect on the correlation functions in the scaling limit of the XYZ chain can be found in Ref. [30].
Note that, as for the QEC model, the effective action is not just the qSG action of Eq. (12), but include corrections, which are irrelevant in the renormalization group sense (see Ref. [30] for explicit forms for these correction terms). However, in contrast to the QEC model, the local spin operator corresponding to the vertex-operator accessible in the XYZ case is e iβφ 2 , which is semi-local. As will be shown later, the expectation value of this vertex operator as well as entanglement properties of the qSG limit of the XYZ chain are more susceptible to corrections to scaling compared to that obtained from the QEC lattice model.

Zero-temperature computations of correlation functions
In this section, we compute various zero-temperature thermodynamic properties of the two different lattice regularizations of the qSG model using DMRG and compare with analytical predictions for the qSG field theory. The DMRG computations were performed using the TeNPy package [65].

QEC model
First, we present the DMRG results for the QEC lattice model. The local Hilbert space on each island was truncated to 9 dimensions: n i = −4, −3, . . . , 3, 4 (similar computations for the free, compactified boson CFT, including the evidence that this truncation of the local Hilbert space is sufficient, are done in Refs. [60,66]). For definiteness, we chose δ = 0.2 and set the overall energy scale by choosing E C 0 = 1 for all the computations of the QEC model.
First, we consider the case E J 0 = 0, when the system is described by a free compactified, boson CFT. The three key properties of this CFT are its central charge (c), the compactification radius R = 1/ √ πK and the plasmon velocity u. They can be obtained from DMRG computations in the following way. The central charge can be verified by computing the scaling of the entanglement entropy (S) with correlation length (ξ) for a partitioning of an infinite system into two semi-infinite halves 4 [67,68,69,70]: The Luttinger parameter of the theory is obtained by computing the algebraic decay of the correlation function [62,63]: Finally, the plasmon velocity u is obtained by computing the zero-temperature ground-state energy (E 0 ) for a finite-system of size L with open (free) boundary conditions and fitting to the Cardy formula [71,72]: where E 0 L is the extensive contribution to the ground state energy. The results are shown in infinite DMRG (left and center panels). We also obtain the ground state energy density E 0 from the infinite system simulations. The right panel of Fig. 2 shows the variation of the ground state energy as a function of system size for free boundary conditions obtained using finite DMRG. By fitting to Eq. (20) and using the obtained value of the central charge, we get the ground-state energy density E 0 and the plasmon velocity u. As shown, both the finite and infinite DMRG results for E 0 match to the third decimal place. The extracted plasmon velocity and the Luttinger parameter values as E J /E C 0 is varied are shown with solid diamonds in Fig. 3. The corresponding perturbative analytical estimates, from Eq. (9), are shown with crosses. While the analytical and the DMRG results approach each other for large E J /E C 0 , for the parameters of interest in this work, the perturbative analytical estimates are clearly insufficient. Now, we consider the case when E J 0 0. First, we compute the expectation value of the local lattice vertex operator e iϕ i and compare with the QFT predictions for the continuum vertex operator e iβφ , given in Eq. (5). Here i is a lattice point within an infinitely large chain. The results obtained with infinite DMRG are shown in Fig. 4. For a given choice E J /E C 0 , which fixes the Luttinger parameter and thus, β [see Eq. (13)], we compute e iϕ i as a function of E J 0 /2E C 0 = M 0 /2. We verify the expected algebraic dependence [see Eqs. (2,3,5)] for different choices of E J /E C 0 (note the log-log scale for the plot). For each choice of E J /E C 0 , the slope is β 2 /(8π − β 2 ). The latter can be used to compute the value of β 2 /8π which is shown in the legend of the plot. In parenthesis, for each curve, the expected value of β 2 /8π is shown obtained by computing the Luttinger parameter for E J 0 /E C 0 = 0 [see Fig. 3]. The agreement is reasonably good, and improves as E J /E C 0 is increased. This improvement is because the corrections to the qSG description of the lattice model, shown in Eq. (10), become more irrelevant in the renormalization group sense.
Note that the above computation only confirms the overall scaling of the expectation value   Here, i is a site within an infinitely large QEC lattice. For these choices of E J /E C 0 , the obtained Luttinger parameters are shown in Fig. 3. The Luttinger parameter determines the qSG coupling β through Eq. (13). The value of β 2 /8π computed using K is shown in parenthesis in the legend for each curve. The algebraic dependence predicted in Eq. (5) is verified, with the slope giving β 2 /(8π − β 2 ). The value of β 2 /8π obtained from the slopes is shown in the legend. As is evident, the agreement is quite good and improves as E J /E C 0 increases. This is because increasing E J /E C 0 makes the irrelevant corrections to the scaling field theory action, given in Eq. (10), less dominant. Then, the qSG description is better suited.
of the local field, which could be inferred purely from dimensional analysis. However, the exact magnitudes of the lattice and continuum vertex operators are equal only up to a proportionality factor, which depends on β or equivalently E J /E C 0 . This proportionality constant arises since the predictions for the qSG field theory assume a certain normalization for the correlation functions of the vertex operators in the conformal limit. 5 We determine this constant of proportionality and use this to compute the two-point correlations of the vertex operators: e iβφ(0) e −iβφ(r) . The corresponding operators to be computed using DMRG are e iϕ i e −iϕ i+r . Using the corresponding form-factors for vertex operators [73,74], we compute the relevant two-point correlation function. The form-factor computation is performed by including contributions up to second breather mass (see Appendix A for more details). Here we only present the comparison between the DMRG and the analytical results in Fig. 5. We chose E J /E C 0 = 1.55 (i.e., β 2 /8π = 0.063) and E J 0 /E C 0 = 0.016 [see also Eq. (13)]. The corresponding soliton mass, determined using Eq. (2) is 0.662. The plasmon velocity determined using the Casimir energy is 1.46, see Fig. 3. The overall field normalization is determined using by computing the expectation value of the one-point expectation value (see Fig. 4). Note that there are no fit parameters in this plot since the mass of the soliton is determined analytically and the field normalization is determined from a different, independent computation. Similar results were obtained for other choices of E J /E C 0 and are not shown for brevity. The infinite DMRG computations are shown with maroon dots, while the form-factor predictions with a solid green line.

XYZ model
In this section, we present DMRG results for the qSG limit of the XYZ chain. We choose the same set of β as in Sec. 4.1 and vary the mass-parameter of the qSG action M 0 [see Eq. (12)]. To make comparison with the QEC lattice model easier, we plot the results with respect to the corresponding QEC mass parameter E J 0 /E C 0 [see Eq. (13)]. For the simulations, we set J x = 1. From Eqs. (2,15,16), the corresponding values of J y , J z are inferred. The results obtained with infinite DMRG are shown in Fig. 6. In the left panel, we plot σ + i ∼ e iβφ/2 as a function of M 0 /2 = E J 0 /2E C 0 on a log-log scale. This expectation value scales as [53] σ + i ∼ e iβφ/2 ∼ The proportionality constant of the second relationship is also conjectured [53], but the exponent of the algebraic dependence is sufficient for our purposes. The extracted values of β 2 /8π are shown for the different data points with the expected values within parenthesis. As an extra check, we compute also the soliton mass (M) as a function of E J 0 /E C 0 . This is directly computed from Eq. (17) by eliminating e iβφ/2 . The results obtained with infinite DMRG are shown on the right panel (empty circles), together with the analytical values of Eq. (16) (crosses). As is evident, the agreement is reasonable, which gives us confidence that we are indeed probing the qSG regime of the XYZ parameter space. . The overall field normalization to relate the lattice operators to the continuum ones is determined by computing the one-point functional (Fig. 4). Note that there are no fit parameters in this plot.   Fig. 4). From the slope, the extracted values of the β 2 /8π [see Eq. (21)] are shown in the legend, together with the expected values in parenthesis. The soliton mass is obtained from the infinite DMRG results for σ + i using Eq. (17). These are shown on the right panel with empty circles for the different choices of β 2 /8π. The corresponding analytical predictions [Eq. (16)] are shown with crosses. Notice that the x-axis values in both panels is much smaller compared to Fig. 4. This is because the XYZ chain is no longer in the scaling regime for the range of E J 0 /E C 0 shown in Fig. 4. This is discussed in Sec. 4.3. Note that the x-axis values for both the panels in Fig. 6 are several orders of magnitude smaller than that of the QEC simulations of Fig. 4. This is because for the values of the QEC simulations, the XYZ chain is no longer in the scaling limit. The reader will notice that the last few data points in the plot of the soliton mass for β 2 /8π = 0.054 (in blue) already start showing deviations occurring due to corrections to scaling. We discuss this next.

Corrections to scaling: QEC vs XYZ
So far, we have computed various zero-temperature thermodynamic quantities for the qSG field theory using DMRG for the QEC and the XYZ lattice regularizations. Now, we discuss how the two regularizations fare with regards to corrections to scaling. Consider a lattice model whose long wavelength behavior (i.e., the scaling limit) is governed by a QFT. Then, there are two corrections to the scaling-limit behavior. First, the lattice Hamiltonian gives rise to various terms in addition to the QFT Hamiltonian, which are irrelevant in the renormalization group sense. The contributions of these terms are small, but nonzero. Second, the QFT operators are approximately represented by the lattice operators. Thus, while computing correlation functions of QFT operators using lattice regularizations, both these corrections contribute to the subleading corrections. For the QEC incarnation of the qSG model, the corrections to the scaling-limit action include the irrelevant term shown in Eq. (10), as well as irrelevant terms formed by higher descendants of the vertex operator e iβφ . A similar set of corrections arise for the effective Hamiltonian of the XYZ chain regularization -some of these terms are explicitly computed in Ref. [30]. However, the two regularizations behave differently when it comes to the definitions of the vertex operators whose correlations are computed. This is because while the spin-operator σ + i of the XYZ chain is approximately equal to e iβφ/2 [see Eq. (17)], the lattice operator e iϕ i of the QEC lattice model, to leading order, is equal to the operator e iβφ up to an overall β-dependent proportionality constant. 6 We observe numerically that larger correlation lengths, ξ XYZ , (equivalently, smaller soliton mass, M), are required to reach the scaling regime of the XYZ spin-chain compared to the QEC lattice. However, this is not always possible in practice due to the following. From Eq. (15), for ξ XYZ → ∞, the parameter l → 0. This is accomplished by choosing |J 2 x − J 2 y | |J 2 x − J 2 z |. However, as β 2 gets close to either 0 or 8π, this becomes increasingly difficult since |J z /J x | → 1 for these choices. This is the other crucial difference between the QEC and the XYZ regularizations. The QEC regularization allows the two qSG parameters: β, M 0 to be independently controlled: the first being controlled by the junction energies (E J ) of the Josephson junctions on the horizontal links, while the latter being controlled by the that (E J 0 ) of the Josephson junctions on the vertical links. In contrast, in the XYZ regularization, M 0 cannot be tuned independent of β.
This difference in corrections to the scaling limits for the two regularizations is shown in Fig. 7. The top left and top right panels show the variation of the expectation values of local fields as a function of the qSG mass-parameter M 0 /2 ∝ E J 0 /2E C 0 [see Eqs. (12,13)]. The QEC lattice results 6 In this discussion, we have concentrated on the qSG vertex operators whose correlation functions are observable by measuring, for instance, current-current correlation functions in a QEC experiment. However, this does not preclude the existence of local operators which have equal or better resilience to corrections to scaling in the XYZ model. Such operators are likely to be superpositions of local qSG operators specifically chosen to be local in the spin language.
for e iϕ i ∼ e iβφ and the XYZ results for σ + i ∼ e iβφ/2 are shown. The qSG field theory predicts a linear-dependence with M 0 for both vertex operators on a log-log scale [see Eqs. (5,17,21)]. As seen from Fig. 7, the QEC vertex operator exhibits this behavior for all choices of M 0 . On the other hand, the corresponding XYZ spin operator does so only deep in the scaling regime when M 0 < 10 −3 for the shown choices of β 2 /8π. One can identify the region where the corrections to scaling are noticeable as the region where |J x − J y | is no longer |J x − J z |. In this region, l is no longer small and thus, the correlation length ξ XYZ is no longer large [Eq. (15)]. Note that the problem is more severe for β 2 /8π being closer to either 0 or 1, which is also apparent from the different curves plotted in Fig. 7 (top right panel). The corresponding correlation lengths are shown in the bottom left and bottom right panels. We do not go beyond E J 0 /2E C 0 = 0.1. This is because for this choice, the correlation is already only a few lattice sites for both models. Thus, beyond this point, the field theory predictions are not expected to describe either lattice model very well. Note that beyond E J 0 /E C 0 > 10 −2 , the correlation length of the XYZ chain actually goes up, while that in the QEC model keeps going down. The increase in ξ XYZ for this range of E J 0 /E C 0 can be understood by noticing that in this case, |J x − J y | ≥ |J x − J z | and the model is no longer in the regime J x > J y ≥ |J z | [see below Eq. (14)]. Recalling that the phase-diagram of the XYZ chain is symmetric under permutation of the coupling constants [31], we can interchange J y and J z and arrive at a corresponding formula for ξ XYZ [Eq. (15)] that is consistent with this behavior. Thus, in this regime, the correspondence between the qSG and the XYZ model presented in Sec. 3.2 is no longer valid. To check that we are not encountering any numerical artifacts of the DMRG simulations, we provide an extra check by performing computations of the entanglement spectrum of the XYZ chain and comparing with analytical predictions also in the regime E J 0 /E C 0 > 10 −2 (see below). Note that there is no such restriction on the parameter space for the QEC model since the two qSG parameters, β, M 0 , can be independently controlled by tuning corresponding lattice parameters E J , E J 0 .

Entanglement spectrum of the qSG model
In the previous section, we have computed various thermodynamic quantities of the qSG model, namely, the one and two-point functions of vertex operators with the QEC regularization and discussed how the scaling regime is reached in comparison to the XYZ chain. Now, we compute the entanglement spectrum of the qSG model.
First, we briefly summarize the generic behavior of the entanglement spectrum for a partitioning of an infinite system into two halves for massive field theories following Ref. [75]. Consider a CFT perturbed by a single primary field Φ (the generalization to multiple perturbations is straightforward). The entanglement spectrum for this massive theory is given by the physical spectrum of a corresponding boundary CFT over a length interval ln ξ, where ξ is the (finite) correlation length [75]. The two spectra are equal up to rescalings and overall shifts and comprise equidistant levels. The boundary CFT has two boundary conditions: free boundary condition at one end, which arises from the entanglement cut and a boundary field, Φ, at the other end (see Fig. 8). The equality of the two spectra is correct up to exponential corrections. This leads to a restriction on the number of low-lying entanglement energy levels which are in correspondence with the boundary CFT spectrum [75]. Note that when the system is critical i.e., Φ = 0 and ξ → ∞, the second F igure 7: Comparison of the expectation values of the vertex operators obtained using infinite DMRG for the QEC lattice model ( e iϕ i ∼ e iβφ , top left) and the XYZ chain ( σ + i ∼ e iβφ/2 , top right) as a function of M 0 /2 ∝ E J 0 /2E C 0 , the mass-parameter of the qSG field theory [see Eqs. (12,13)]. The field theory computations for both the quantities predict a linear dependence (see Figs. 4,6). For very small M 0 /2, where the correlation-length is large for both QEC and the XYZ lattice models and the continuum qSG description is valid for both models. As E J 0 /E C 0 increases, the correlation length diminishes and the corrections to scaling becomes important. However, as seen from the QEC plot, the expectation value of the vertex operator continues to follow the straight line, while that from the XYZ chain deviates from the field theory predictions at E J 0 /E C 0 ∼ 10 −3 . This difficulty of reaching the scaling regime for the XYZ chain occurs in the region where |J x − J y | is no longer |J x − J z |. Then l no longer tends to zero and thus, the correlation length ξ XYZ is no longer large [see Eq. (15)]. We do not go beyond E J 0 /2E C 0 = 0.1 since the correlation length for both the XYZ and the QEC models is only a few lattice sites and it is not meaningful to apply the qSG field theory predictions beyond this point. For reference, the corresponding correlation lengths are shown in the bottom left and bottom right panels. Note that the correlation length for the XYZ chain actually goes up for E J 0 /E C 0 > 10 −2 . This is because for this choice of parameters, the XYZ chain is no longer in the scaling limit and increasing the aforementioned ratio no longer corresponds to increasing the mass gap of the model (see main text for more details). boundary condition (that on the left end in Fig. 8) is inherited from the original model whose entanglement spectrum is being computed [60,76].
The aforementioned relationship between the two spectra holds independently of whether the perturbation, Φ, is integrable or not. But, the qSG model is an integrable deformation of the free, compactified boson CFT and thus, the question arises as to whether one can use integrability to glean additional information about the qSG entanglement spectrum. To that end, we can use the fact that the qSG model arises as the scaling limit of the quantum XYZ chain (Sec. 3.2). The latter spin-chain or its "classical" version, the eight-vertex model [31], falls within the category of integrable lattice models which exhibit equidistant levels for the entire entanglement spectra -notable other examples include the transverse-field Ising and the XXZ models [77,78,79]. We compute the entanglement spectrum of the XYZ chain both analytically and numerically in Sec. 5.1. In particular, we show that the level-spacing, denoted by ε XYZ , goes as 1/ ln ξ XYZ , where ξ XYZ is the correlation length of the XYZ model [see Eq. (15)] as long as the system size is much larger than ξ XYZ . This behavior is consistent with what is predicted in Ref. [75].
After computing the qSG entanglement spectrum using the XYZ chain, we compute the same for the QEC lattice model using DMRG. The primary motivation for this computation is to investigate to what extent the low-lying entanglement spectrum of the XYZ chain is a universal feature of the qSG model. As will be shown in Sec. 5.2, the spectra computed using the XYZ and the QEC lattices have identical degeneracies. Furthermore, the level spacing of the entanglement spectrum computed using the QEC lattice model, denoted by ε QEC also goes as 1/ ln ξ QEC , where ξ QEC is the correlation length fo the QEC lattice model. Without fine-tuning of the microscopic models, there is no reason why ε XYZ would be equal to ε QEC . However, both quantities scale inversely with the logarithm of the respective correlation lengths. At the same time, both ξ XYZ and ξ QEC depend on the qSG mass-parameter, M 0 , as M −1/(2−β 2 /4π) 0 . Thus, from purely dimensional considerations, we can conclude that ε XYZ and ε QEC are linearly dependent on each other: where a 0 (β), a 1 (β) are (possibly non-universal) functions that depend on the parameters of the two lattice models. Here, we have also explicitly indicated the dependence of the entanglement level spacings on the qSG parameters: β, M 0 . We verify this linear dependence in Sec. 5.2. The secondary motivation for this computation is to demonstrate that the QEC lattice model continues to provide meaningful prediction for the qSG entanglement spectrum even when the XYZ chain is no longer in the scaling limit (see discussion of Sec. 4.3).

XYZ model
Now, we compute the qSG entanglement spectrum using the XYZ chain. The XYZ Hamiltonian [see Eq. (14)] can be related to the transfer-matrix of the classical eight-vertex model. To establish this relation, it is useful to consider the principal regime for the two models [31]. Denote the XYZ couplings in the principal regime by J p α , where α = x, y, z. Then, the principal regime is given by The couplings of the XYZ chain can be related to the two parameters of the classical eight-vertex model, denoted by Γ, ∆. In the principal regime, they are given by As a result, in the principal regime, for the eight-vertex model, Since we are interested in the entanglement properties of the XYZ chain, we need two further parameters, λ, k, given by [31] 2 where sn is the Jacobi sine function [31]. Here, 0 ≤ k ≤ 1 and 0 ≤ λ ≤ I(k ), where k = √ 1 − k 2 and I(k) is the complete elliptic integral of the first kind with modulus k. The entanglement spectrum for the XYZ chain can be related to the spectrum of the CTM of the eight-vertex model [32,33]. The entanglement spectrum comprises equidistant levels, with the level spacing given by ε XYZ , given by [31] Now, we compute the scaling of the level-spacing of the entanglement spectrum, ε XYZ , as the mass-parameter of the sine-Gordon action, M 0 of Eq. (1), is taken to zero. This corresponds to taking the limit J y /J x → 1 − in the XYZ chain. For the purposes of the calculation, we set J x = 1 and consider the case when |J y /J z | ≥ 1 (the other case can be analyzed similarly). Our goal is to compute the scaling of the level-spacing of the entanglement spectrum, ε XYZ , given in Eq. (27) as J y /J x → 1 − . To that end, we first define the couplings in the principal regime: Thus, the limit J y /J x → 1 − is equivalent to J p x /J p z → −1 + , which in turn implies ∆ p → −1 − [see Eq. (24)]. From Eq. (26), this implies k → 1 − . We will also use the fact that −isn(iλ, k) → tan λ, I(k) → − 1 2 ln 1 − k 8 (29) as k → 1 − . To quantify deviations from the critical point, we define two small parameters: Then, from Eq. (26), we get where we have kept terms up to O(x). Next, we use l 2 2x sin 2 (β 2 /8) as k → 1 − to get which is commensurate with the general statement that the entanglement spectrum gap closes as 1/ ln ξ XYZ to leading order [75,79], see Eq. (15). We check the leading order term by analyzing the dependence of ε XYZ as a function of l. This is shown in Fig. 10, which confirms the predicted linear dependence with 1/ ln(l/4) with a slope that is close to −π 2 (1 − β 2 /8π).

QEC model
In this section, we compute the qSG entanglement spectrum using the QEC lattice model. As discussed earlier, the low-lying part of the spectrum should be a characteristic of the qSG field theory and thus, should be comparable to the results obtained from the XYZ chain (see in Fig. 9). The results are shown in Fig. 11 for β 2 /8π 0.063 (similar results were obtained for other choices and are not shown for brevity). As seen from the left panel, the low-lying entanglement spectrum exhibits the expected equidistant level structure, with the same degeneracy structure given in Fig. 9 (a). Despite the overall degeneracy structure being consistent, it is clear that data quality for the QEC model is worse compared to the XYZ chain. One of the reasons for this is that the XYZ chain, unlike the QEC model, exhibits the equidistant structure for the entire entanglement spectrum due to its relationship to the eight-vertex model. We are not aware of any such deep connections for the QEC model. We believe the worse data quality is also partially due to the non-universal lattice effects which affect the two models differently. Finally, at a more pragmatic level, the large local Hilbert space of the QEC model makes the computations much more resource-consuming compared to the XYZ chain. This restricts the size of the bonddimensions that are accessible for a moderate-scale simulation effort pursued in this work. In the center panel, we show the level-spacings for different choices of E J 0 /E C 0 as obtained from the QEC and the XYZ models for β 2 /8π = 0.063. As argued earlier [see discussion before Eq. (22)], the two level-spacings are not equal to each other, but are linearly related. The verification of this linear dependence is shown in the right panel of Fig. 11.   Next, in Fig. 12 (left panel), we verify the linear dependence of ε QEC on ε XYZ for different choices of β 2 /8π. In the right panel, the corresponding parameters of the linear dependencies, given in Eq. (22), are shown. At this point, we do not have a deep understanding of the a 0 , a 1 and their functional dependency on β. It is plausible that these are non-universal functions of β which depend on the details of the QEC and XYZ lattice models, but we leave a detailed investigation for a future work.

Corrections to scaling in the entanglement spectrum: QEC vs XYZ
Here, we demonstrate that the corrections to scaling, which are noticeable for E J 0 /E C 0 > 10 −3 and lead to incorrect dependence of the qSG vertex operator e iβφ/2 in the XYZ chain (see Sec. 4.3), also causes the qSG entanglement spectrum to be incorrectly inferred from the XYZ results. This is in contrast to the QEC model, which continues to provide meaningful physical predictions for the entanglement spectrum for these choices of E J 0 /E C 0 . The results are shown in Fig. 13 for β 2 /8π 0.063 (similar results were obtained for other values and are not shown for brevity). For E J 0 /2E C 0 < 10 −3 , both XYZ and the QEC models provide meaningful results for the qSG entanglement spectrum and the level-spacings, ε QEC and ε XYZ , are linearly dependent on each other (see Fig. 12). However, upon further increase of E J 0 /E C 0 , ε XYZ grows much faster violating the linear dependence. The point of departure of the linear dependence coincides precisely with the departure of the linear dependence of the vertex operator e iβφ/2 ∼ σ + in Fig. 7. Beyond E J 0 /2E C 0 ∼ 10 −2 , ε XYZ actually decreases. Clearly, for these parameters, ε XYZ cannot correspond to the levelspacing of the qSG entanglement spectrum. This is because increasing E J 0 /E C 0 increases the mass-gap of the qSG model, which would increase the entanglement-spectrum level spacing of the   Fig. 7 (bottom right panel), this region of decreasing ε XYZ corresponds precisely to the region where the correlation length of the XYZ chain, ξ XYZ , increases. We again emphasize that the physics of the XYZ model is perfectly consistent -increasing correlation length should be associated with a decreasing entanglement level-spacing. An independent check of this is provided by the CTM calculations of the entanglement spectrum (see Sec. 5.1). As seen from Fig. 13, the agreement between the DMRG and the CTM calculations is excellent for all choices of E J 0 /E C 0 . Finally, the results obtained using the QEC lattice exhibit the expected behavior for all choices of E J 0 /E C 0 . Before concluding this section, we point out that the qSG entanglement spectrum could, in principle, be inferred using the spectrum of the boundary sine-Gordon model [42,80,81] after taking the strength of the bulk perturbation of the latter model to zero. Potentially, this could also be a way to investigate the fate of the boundary bound-states of the boundary sine-Gordon model in the massless bulk limit and we hope to return to this problem in the future. We did check that the degeneracies of the entanglement spectrum of the qSG do indeed match that of the free, compactified boson CFT with Dirichlet boundary conditions [60]. The latter boundary condition can be viewed as an extreme case when the strength of the boundary potential is taken to infinity.

Summary and outlook
To summarize, we numerically analyzed with DMRG a faithful, analog, quantum simulator built with QEC elements for the qSG model in 1+1 space-time dimensions. The QEC model provides a lattice-regularization of the qSG model using local interactions that can be physically  Figure 13: Level-spacing of the qSG entanglement spectrum as computed with the QEC and the XYZ lattice models. We chose β 2 /8π 0.063 (similar results were obtained for other choices and are not shown for brevity). For both models, we show the DMRG results. Furthermore, for the XYZ chain, we show the analytic CTM predictions (see Sec. 5.1) as well. For E J 0 /E C 0 < 10 −3 , both the models give physically meaningful predictions, with ε QEC being linearly dependent on ε XYZ . However, further increase of E J 0 /E C 0 causes ε XYZ to increase much faster compared to ε QEC . Beyond E J 0 /E C 0 > 10 −2 , ε XYZ goes down, which is clearly incompatible with the expectation that the entanglement level-spacing increases with increase of the qSG mass-parameter. The results for the QEC model, however, continue to provide physically meaningful predictions, increasing steadily with E J 0 /E C 0 . The perfect overlap of the CTM and the DMRG results for the XYZ chain show that what we are observing is not a numerical artifact of DMRG; rather, it is the qSG-XYZ correspondence that is no longer valid for E J 0 /E C 0 > 10 −3 (see Secs. 4.3, 5.3 for more discussion).
realized by a straightforward generalization of current experimental works. By computing various zero-temperature thermodynamic properties of the QEC model with DMRG and comparing with the qSG field theory predictions, we numerically demonstrate that the QEC lattice indeed realizes the qSG model. Furthermore, we show that in contrast to the integrable XYZ-chain regularization, the QEC lattice model is less susceptible to corrections to scaling. In contrast to the XYZ chain, where the spin-operator σ + corresponds to the qSG vertex operator e iβφ/2 , the QEC model directly starts with lattice versions of the operators e iβφ . Furthermore, we computed the entanglement spectrum of the qSG model using both the XYZ and the QEC models and showed that the lowlying entanglement energy levels exhibit the same set of degeneracies. We provided a scaling argument to show that the level-spacings of the low-lying entanglement spectrum for the two models are linearly related to each other and verified this claim with numerical results. Finally, we also showed that in the XYZ chain, the same corrections to scaling that plague the correlation functions of the vertex operators also cause the model to predict unphysical values of the qSG entanglement spectrum. The latter problem is also remedied by the QEC lattice model. The current work gives rise to many, new, potentially-fruitful research directions and we discuss some of them below. First, concerning the qSG model, an experimentally-realizable, numerically-tractable lattice model potentially opens the door to the investigation of several open problems -examples include many finite temperature properties of the correlation functions for the qSG model [82,83,84]. At the same time, this work opens the possibility of experimentally exploring non-equilibrium phenomena in the qSG model, which, in the recent years, have received enormous interest [85,86,87,88,89,90,91,92,93]. Second, our approach to faithfully simulate QFTs with QECs can be used to investigate multi-field models. For the latter, it is crucial that the underlying lattice degrees of freedom faithfully give rise to the continuum ones, without resorting to mathematical manipulations like bosonization. This is because properties as fundamental to the QFT as integrability can be difficult to relate in the fermionic and bosonic counterpartse.g., the quantum double sine-Gordon model, which can be faithfully realized with QECs [13]. We aim to analyze the latter model with QECs in the near future. Third, the analog QEC simulator for the qSG model analyzed in this work can be readily generalized to include integrability-breaking perturbations [13]. After all, the eventual goal is to simulate interacting QFTs to answer questions which are intractable with analytical methods. This is possible with QEC lattices, which can be used to simulate the massive Schwinger model [94] or a two-frequency generalization of the qSG model [95,96]. In fact, the basic primitives for realizing the two-frequency generalization have already been experimentally demonstrated [97,98]. Fourth, analog QEC simulation of QFTs can also be implemented in higher dimensions. In particular, QEC arrays in 2+1 space-time dimensions have already been built and analyzed experimentally in the context of realizing interacting bosonic models, even in hyperbolic space [99]. Given this recent experimental progress, we are optimistic of the use of QECs for investigation of 2+1D QFTs in the near future. In this appendix we derive analytic results for the two-point correlation function e iϕ i e −iϕ i+r shown in Fig. 5. More precisely we calculate the static, zero-temperature two-point function e iβφ(0) e −iβφ(r) via a form-factor expansion [73] directly in the continuum model (1). The basic idea is to insert a resolution of the identity between the operators, where the sum runs over all possible intermediate states. Since the spectrum of the qSG model is exactly known, these intermediate states can be classified by their particle contents (solitons, antisolitons and breathers of type n) and the respective momenta of the particles. Since the masses of the particles in the intermediate state will lead to an exponential decay at large distances, the leading behavior will be governed by the lightest particles. Here we consider the vacuum state |0 , single breathers of type 1 and 2, |θ 1,2 , and two 1-breather states |θ 1 , θ 2 1,1 respectively, where we parametize the momenta of the intermediate n th breather via their rapidities, P = m n u sinh θ. Thus we evaluate e iβφ(0) e −iβφ(r) = 0|e iβφ |0 2 + dθ 2π 0|e iβφ |θ 1 2 e −i m 1 u r sinh θ + dθ 2π 0|e iβφ |θ 2 2 e −i m 2 u r sinh θ + 1 2 dθ 1 dθ 2 (2π) 2 0|e iβφ |θ 1 , θ 2 1,1 2 e −i m 1 u r i sinh θ i + . . .

Acknowledgments
where the factor 1 2 in the last term avoids double counting, and the dots represent terms with heavier intermediate states. For the parameters of Fig. 5, β 2 /8π = 0.063, the next terms would be the single 3-breather state (mass m 3 ), the 1-breather-2-breather pair (mass m 1 + m 2 ) and the three 1-breather state (mass 3m 1 ). The form factors (matrix elements) appearing in the expansion are known exactly [74]. With this a straightforward calculation yields the result e iβφ(0) e −iβφ(r) = G 2 β 1 +