# Quantum electronic circuit simulation of generalized Sine-Gordon models

Ananda Roy, Hubert Saleur

<figure><img src="image_1.png" /></figure>

To cite this version:

Ananda Roy, Hubert Saleur. Quantum electronic circuit simulation of generalized Sine-Gordon models. Phys.Rev.B, 2019, 100 (15), pp.155425. 10.1103/PhysRevB.100.155425. hal-02066928

HAL Id: hal-02066928

https://hal.science/hal-02066928

Submitted on 12 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Quantum electronic circuit simulation of generalized sine-Gordon models

Ananda Roy<sup>1,2,*</sup> and Hubert Saleur<sup>1,3</sup>

<sup>1</sup>Institut de Physique Théorique, Paris Saclay University, CEA, CNRS, F-91191 Gif-sur-Yvette, France

<sup>2</sup>Department of Physics, T42, Technische Universität München, 85748 Garching, Germany

<sup>3</sup>Department of Physics and Astronomy, University of Southern California, Los Angeles, California 90089-0484, USA

<figure><img src="image_3.png" /></figure>

(Received 6 March 2019; revised manuscript received 7 October 2019; published 23 October 2019)

Investigation of strongly interacting, nonlinear quantum field theories (QFTs) remains one of the outstanding challenges of modern physics. Here, we describe analog quantum simulators for nonlinear QFTs using mesoscopic superconducting circuit lattices. Using the Josephson effect as the source of nonlinear interaction, we investigate generalizations of the quantum sine-Gordon model. In particular, we consider a two-field generalization, the double sine-Gordon model. In contrast to the sine-Gordon model, this model can be *purely quantum integrable*, when it does not admit a semiclassical description—a property that is generic to many multifield QFTs. The primary goal of this work is to investigate different thermodynamic properties of the double sine-Gordon model and propose experiments that can capture its subtle quantum integrability. First, we analytically compute the mass spectrum and the ground-state energy in the presence of an external “magnetic” field using Bethe ansatz and conformal perturbation theory. Second, we calculate the thermodynamic Bethe ansatz equations for the model and analyze its finite temperature properties. Third, we propose experiments to verify the theoretical predictions.

DOI: 10.1103/PhysRevB.100.155425

I. INTRODUCTION

The longstanding goal of quantum field theory (QFT) is to predict the masses of the excitations and their scattering cross sections in terms of the parameters of the theory and to characterize the different phases and the phase-transition points. While remarkable progress has been achieved in analyzing strongly coupled QFTs using numerical or effective field theory methods, many quantities of interest remain elusive, either due to the computational limitations or due to the lack of an effective field theory. Quantum simulation, both analog and digital, takes a different approach to solving these aforementioned problems, where one quantum system is tailored to simulate another in a controlled manner [1-9].

Here, we describe analog quantum simulators for some QFTs with superconducting quantum electronic circuit (QEC) lattices [10,11]. Specifically, we are interested in obtaining bosonic QFTs directly, and not as consequences of mathematical manipulations (bosonization) of an underlying fermionic or spin system [12-17]. Even in 1 + 1 space-time dimensions, this is crucial to distill the true behavior of a bosonic system from a bosonized fermionic system [18,19], despite the well-known fermion-boson correspondences [20,21]. Of course, this distinction is even more important in higher dimensions.

Attempts in this direction have been considered previously in the well-established Bose-Hubbard model paradigm [22,23]. Another scheme is to use the Josephson effect to give rise naturally to the cosine potential [24] of the sine-Gordon (SG) model (see more below). Our construction generalizes the latter strategy. The role of the bosonic field at a point in

space-time is played by the time integral of the voltage at that point [25] (so the underlying degrees of freedom of the QEC lattices are faithful directly to the bosonic description). The building blocks of the QEC lattices are superconducting self- and mutual inductances, capacitances, and Josephson junctions. The current-voltage constitutive relations of these elements, together with Kirchhoff's laws for the circuits, give rise to the nonlinear field equations of the QFTs.

The playground for bosonic nonlinear QFTs is wide, even in 1 + 1 dimensions, and includes typically cousins of the SG theory involving multiple component fields, together with some cosine interaction. The case of two bosons already escapes our complete understanding. In contrast to the SG case, the physics of the various “double sine-Gordon models” (see definitions below) cannot always be inferred from some classical limit. For instance, fixed points are known to exist, which are truly quantum in nature and do not admit a semiclassical description [26]. In fact, one of the simplest new aspects to formulate and study in multicomponent bosonic theories is integrability. Integrability leads to a factorized multiparticle scattering matrix [27,28]. This, in turn, allows analytic computation of thermodynamic and transport properties of these QFTs [29,30]. The ability to do such computations is crucial to understand nonperturbative aspects of strongly interacting systems—for instance, charge fractionalization [31-33] and its effect on full-counting statistics [34], both topics of high current interest.

It is well known that the bulk and boundary SG models are classically integrable [35,36]. This was largely used historically to establish their quantum integrability as well [28,37,38]. While the naive intuition that quantum fluctuations might well destroy classical integrability—as is known to be the case, for instance, for many sigma models [39]—

<footnote>* ananda.roy@tum.de</footnote>
is usually correct, it is now understood that the SG model remains integrable in the quantum regime due to the existence of subtle quantum group symmetries [40]. The presence of these symmetries in 1+1 dimensions relies heavily on the existence of conformally invariant UV fixed points, of which integrable QFTs can be considered as deformations [41–45].

The simplest generalization of the SG model in this regard is probably the double sine-Gordon (dSG) model, whose Euclidean action is

$$S_{\text{dSG}} = \int d^2x \frac{1}{2} \sum_{i=1,2} (\partial_\mu \phi_i)^2 + \frac{2M_0}{\pi} \prod_{i=1,2} \cos(\alpha_i \phi_i), \quad (1)$$

where $\phi_1, \phi_2$ are two bosonic fields, $M_0$ is the interaction strength, and $\alpha_1, \alpha_2$ are the coupling constants. This model and its close cousins involving a boundary interaction appear in a variety of guises in condensed-matter physics—for instance, in transport experiments involving one-dimensional particles with charge and spin, or in quantum Brownian motion on two-dimensional lattices [46]. Unfortunately, little remains known about the possibility to solve the dSG model exactly, and to infer from such a solution properties of physical interest. This is because, remarkably, as soon as more than one bosonic degree of freedom is involved, quantum and classical integrability often part ways. In fact, it is known that the dSG model is classically integrable only for $\alpha_1 = \alpha_2$, when it reduces to two decoupled SG models [47]. However, the same fluctuations that can destroy classical integrability in some cases can also give rise to quantum integrability. Purely quantum integrable manifolds, if they exist, should thus arise for values of the coupling constants that cannot approach the origin in the $(\alpha_1, \alpha_2)$ plane. The dSG model is probably only quantum integrable on several $(\alpha_1, \alpha_2)$ manifolds [15,18,48–51]. While some doubts remain on the exact nature of this statement, it is strongly believed in the community that the manifold, $\alpha_1^2 + \alpha_2^2 = 4\pi/\hbar$, indeed is quantum integrable (in the following we will set $\hbar = 1$). This is, of course, a remarkable statement: we are facing a situation where a classical soliton wave packet, which solves the classical field equations, gets scrambled, but its quantized counterpart, when the couplings reach some magical values, propagates undistorted and scatters only with phase shifts! One of our goals is to propose a setup where such behavior—which, we believe, is the norm rather than the exception in many multifield theories—might be observed experimentally.

It is important to stress that the calculations presented in this paper are done exactly for the model in Eq. (1). It should be contrasted with those done by extrapolating the results of the more general Fateev model [48], of which the dSG is a special case. However, the Fateev model involves, on top of the two (compact) bosonic fields $\phi_i$, an extra noncompact bosonic field. The presence of this other field makes many calculations of the Fateev model—in particular, the thermodynamic Bethe ansatz (TBA)—considerably easier. It is expected that the final results can then be continued to the limit where this extra boson decouples, but this assumption is dangerous in view of the singularity of this limit.

## II. DOUBLE SINE-GORDON MODEL AT ZERO TEMPERATURE

The solitons of the quantum dSG model cannot be approached with semiclassical methods, and are not fully understood. It is known that they are topological excitations of the fields $\phi_1, \phi_2$, where $\phi_{1,2} = (\alpha_1\phi_1 \pm \alpha_2\phi_2)/2\sqrt{\pi}$. They carry a pair of quantum numbers corresponding to the fields $\phi_1, \phi_2$, just as an electron carries electric charge and spin quantum numbers. These solitons scatter with a factorized scattering matrix. More remarkably, because of the underlying (pair of) quantum group symmetries [49], the two quantum numbers of each soliton scatter independently of each other. The factorized scattering matrix is expected to be given by [48] $S = S_{p_1} \otimes S_{p_2}$, where $S_{p_i}$, $i = 1, 2$ is the SG scattering matrix. The corresponding SG couplings $\beta_i$'s are related to $p_i$'s by $\beta_i^2/8\pi = p_i/(p_i + 1)$, $i = 1, 2$ [28,45], where $p_{1,2} = \alpha_{1,2}^2/2\pi$ [48].

The masses of the physical excitations in terms of the parameters of the action, together with a verification of the scattering matrix, are obtained by calculating the ground state energy in the presence of “magnetic fields” $h_{1,2}$ coupling to the conserved charges $Q_{1,2} = \int dx \partial_x \phi_{1,2}/\sqrt{\pi}$. Consider the limit $h_{1,2} \gg M_0$ and tune $h_2$ such that $\langle Q_2 \rangle = 0$. Standard Bethe ansatz calculations using the Wiener-Hopf technique [48,52] show that, on the integrable manifold, the ground-state energy gets a logarithmic correction in $h_1$ (to be contrasted with the purely polynomial corrections of the Fateev model [48]), in addition to its usual quadratic dependence:

$$\mathcal{E}(h_1) = -\frac{h_1^2 p_1 (2 - p_1)}{2\pi} + \frac{m_s^2}{2\pi} \sin^2\left(\frac{\pi p_1}{2}\right) \ln h_1. \quad (2)$$

Here $m_s$ is the mass of the $\phi_1$ soliton [53]. It is obtained in terms of the action parameters using conformal perturbation theory (CPT) [54]. We find that $m_s = M_0/\sin(\pi p_1/2)$, which matches the result of Ref. [48]. The calculation is described in detail in Appendix A.

## III. DOUBLE SINE-GORDON MODEL AT NONZERO TEMPERATURE

The thermodynamics of QFTs in an infinite volume can be inferred from its scattering matrix [55] using the TBA technique [29]. Using the factorized scattering matrix given above, we now compute the free energy of the dSG model at a finite temperature $1/R$. The TBA given below is essential to correctly predict the dSG free energy and cannot be directly inferred from the TBA of the Fateev model (the free energy extrapolated from the Fateev model will appear to be infinite). Without loss of generality, we choose $p_1 < 1$, i.e., the first SG sector to be attractive, and $p_2 > 1$, i.e., the second, repulsive. In addition to the solitons, the spectrum of the theory also includes $n < 1/p_1$ particles which are neutral with respect to the first quantum number, but still carry the second quantum number. They are the bound states of the fundamental “double” solitons and correspond to the ordinary breather poles of the attractive SG scattering amplitude of the first sector [28]. The TBA analysis is, in general, highly nontrivial, the scattering in each sector being nondiagonal. We consider the simpler case $p_1 = 1/(n+1)$ when the first sector scattering
<figure><img src="image_4.png" /><figcaption>FIG. 1. TBA diagram for the dSG model. The solid circles denote the physical massive particles: soliton (s), antisoliton (a), and breathers ($b_j$, $j = 1, ..., n$). The empty circles denote pseudoparticles ($j = 1, ..., n + 2$) needed to diagonalize the dSG scattering matrix. The pseudoparticle 1 has a cross on it, which indicates that this particle has a mass term in its TBA equation. The connectivity of the diagram encodes which particles show up in the TBA equation of a given particle. The arrows on the links encode the sign of the term on the right-hand side (RHS) of the TBA equation. For instance, for a link connecting particles p, q, if there is an arrow incident on q and none on p, then the RHS of the TBA equation for p has a minus sign in front of the term involving q, while the term on the RHS for q involving p has a plus sign.</figcaption></figure>

is diagonal. The diagonalization problem of the second sector is done by the Takahashi-Suzuki classification of the solution of the Bethe equations [56–58] and the algebraic Bethe ansatz [59,60]. The detailed derivation of the TBA equations is given in Appendix B. The resulting TBA equations have a universal form for all n:

$$
\begin{align*}
\epsilon_{\gamma}(\theta) &= \hat{\phi}_n \star \mathcal{L}_{b_n}, & \gamma &= s, a, \\
\epsilon_{b_n}(\theta) &= \hat{\phi}_n \star (\mathcal{L}_s + \mathcal{L}_a + \mathcal{L}_{b_{n-1}}), & \\
\epsilon_{b_j}(\theta) &= \hat{\phi}_n \star (\mathcal{L}_{b_{j+1}} + \mathcal{L}_{b_{j-1}}), & j &= n-1, \dots, 3, \\
\epsilon_{b_2}(\theta) &= \hat{\phi}_n \star (\mathcal{L}_{b_3} + \mathcal{L}_{b_1} - \mathcal{L}_1), & \\
\epsilon_{b_1}(\theta) &= \hat{\phi}_n \star (\mathcal{L}_{b_2} - \mathcal{L}_2), & \\
\epsilon_1(\theta) &= -m_{b_1} R \cosh \theta + \hat{\phi}_n \star (\mathcal{L}_{b_2} - \mathcal{L}_2), & \\
\epsilon_2(\theta) &= -\hat{\phi}_n \star (L_{b_1} - L_1 + L_3), & \\
\epsilon_j(\theta) &= -\hat{\phi}_n \star (L_{j+1} + L_{j-1}), & j &= n-1, \dots, 3, \\
\epsilon_n(\theta) &= -\hat{\phi}_n \star (L_{n-1} + L_{n+1} + L_{n+2}), & \\
\epsilon_k(\theta) &= -\hat{\phi}_n \star L_n, & k &= n+1, n+2. \tag{3}
\end{align*}
$$

The TBA kernel is $\hat{\phi}_n(\theta) = (n+1)/\cosh[(n+1)\theta]$, while $\mathcal{L}(\theta) = \ln[1 + e^{\epsilon(\theta)}]$ and $L(\theta) = \ln[1 + e^{-\epsilon(\theta)}]$. Finally, $m_{b_1} = 2m_s \sin \pi p_1 / 2$ is the mass of the first breather [28] and $a \star b = \int a(\theta - \theta')b(\theta')d\theta'/2\pi$. The TBA diagram is given in Fig. 1.

We check our TBA equations by analytic computation of the central charge in the conformal (UV) limit. This yields the expected result of $c_{UV} = 2$ for the two bosonic fields (see Appendix B for details of the computation). Next, we compute the free energy of the system upon deviation from the conformal limit. This is given in terms of the effective central charge, $c_{eff} = -6R^2f/\pi$, where f is the free energy per unit length. On the integrable manifold, the perturbation has dimension $(\alpha_1^2 + \alpha_2^2)/(4\pi) = 1$. This leads to logarithmic corrections to $c_{eff}$ [61,62] upon deviations from the conformal limit, in addition to the polynomial corrections in powers of $(m_s R)^2$ [54]. We calculate these corrections by numerically solving the TBA equations. The logarithmic correction for the

<figure><img src="image_5.png" /><figcaption>FIG. 2. (a)–(d) Essential primitives of QEC lattices [(a) an inductance, L, (b) a capacitance, C, (c) a mutual inductance, M, and (d) a Josephson junction (junction energy $E_J$ and junction capacitance $C_J$)]. The magnetic fluxes going through the inductance (L), the mutual inductance (M), and the Josephson junction are denoted by $\phi$, $\phi_1$, $\phi_2$, and $\phi$ respectively. The charges on the capacitor plates of the capacitances C and $C_J$ are denoted by Q. (e) One unit cell of the QEC for the dSG model. The top and bottom horizontal lines denote the electrical ground. The Josephson junctions (junction energy $E_{J,0}$ and junction capacitance $C_0$) on the vertical links give rise to the nonlinear interaction. The Josephson junctions (junction energy $E_J$ and junction capacitance $C_J$) give rise to a large impedance of the array, without any nonlinear effects. Finally, the mutual inductance M and the capacitance C provide the necessary inductive and capacitive coupling between the upper and lower parts of the array. The bosonic fields, $\phi_1$, $\phi_2$, are the node fluxes at the points indicated.</figcaption></figure>

different values of n is then verified using CPT. We find the logarithmic correction to $c_{eff}$ to be

$$ (\Delta c_{eff})_{\ln} = \frac{6(m_s R)^2}{\pi^2} \sin^2 \left[ \frac{\pi}{2(n+1)} \right] \ln(m_s R). \quad (4) $$

More details on the TBA and the explicit forms of the different polynomial corrections to $c_{eff}$ for $n=2, 3$ are given in Appendix B. Note that this logarithmic correction is very similar to those obtained for the sine-Gordon model for suitably chosen coupling constants [61,62] and should be contrasted with a very different type of logarithmic correction that arises in the Fateev or sinh-Gordon models [48,63].

## IV. DOUBLE SINE-GORDON MODEL WITH QUANTUM CIRCUITS

The predictions for the various thermodynamic quantities computed above can, in principle, be measured in an experimental setup. In particular, experiments should be able to capture the subtle quantum vs classical integrability of the dSG model. To that end, we provide a QEC for the dSG model. In terms of the rotated fields $\phi_{1,2}$, the dSG action can be viewed as two SG actions coupled by $B \int d^2x \partial_\mu \phi_1 \partial_\mu \phi_2$, where $B = \pi/\alpha_1^2 - \pi/\alpha_2^2$. Recall that in QECs, the bosonic field is the time integral of the voltage at a point. Thus, the first coupling term, $\partial_t \phi_1 \partial_t \phi_2$, is a voltage-voltage coupling, realized in QECs with a capacitance (C). The second coupling term, $\partial_x \phi_1 \partial_x \phi_2$, is a current-current coupling, realized by a mutual inductance (M). All that is left is to realize two separate identical SG models. The QECs for the two SG models are two one-dimensional arrays [see Fig. 2(e)] with Josephson junctions on horizontal links (junction energy $E_J$ and junction
capacitance $C_J$) together with Josephson junctions on vertical links (junction energy $E_{J,0}$ and junction capacitance $C_0$). We work in the regime when $E_J \gg E_{C_J}$ and $E_{C_0} \gg E_{C_J}$, where $E_{C_J} = (2e)^2/2C_J$ and $E_{C_0} = (2e)^2/2C_0$. The array impedance is $Z = R_Q\sqrt{2E_{C_0}/E_J}/2\pi$ [64–67], where $R_Q(\approx 6.5 \text{ k}\Omega)$ is the resistance quantum (in units of Cooper pairs). Only the kinetic inductance of the Josephson junctions on the horizontal links are used to increase the impedance of the array, with the phase-slip amplitudes exponentially small, $\sim e^{-\sqrt{E_J/E_{C_J}}}$ [68]. The Luttinger parameter is $K = 2Z/R_Q$. Capacitively and inductively coupling these two SG models results in the QEC for the dSG model. Choosing $M = C$, the correspondences between the circuit components and the dSG action are $\alpha_1^2 = \pi K/2$, $\alpha_2^2 = \pi/(2K + 4\pi C)$ and $M_0 = \pi E_{J,0}$. A more detailed analysis of the circuit is given in Appendix C.

Next, we describe various experimental measurements of the dSG model possible in a QEC setup. Estimation of $\alpha_1, \alpha_2$ can be done by biasing with two current sources (the “magnetic fields” $h_1, h_2$ considered earlier) and measuring the transmission of an input field through the dSG lattice. This provides an estimate of the average current, $\langle \partial_x \varphi_{1,2} \rangle$, flowing through the horizontal links in the circuit. Comparison to CPT calculations, which give $\langle \partial_x \varphi_{1,2} \rangle = [(\alpha_1 \pm \alpha_2)h_1 + (\alpha_1 \mp \alpha_2)h_2]/4\pi^{3/2}$ to leading order, yields estimates of $\alpha_1, \alpha_2$. Classical integrability can be captured by scattering measurements of a classical sine-Gordon soliton wave packet [36] for the either $\varphi_1$ or $\varphi_2$ propagating through the array. Due to the presence of the mutual inductance and capacitance coupling to the upper and lower part of the array (Fig. 2), the wave packet will distort. However, on the quantum integrable manifold, a quantum soliton would propagate undistorted. A signature of the factorized scattering of the quantum soliton can be obtained by measuring the specific heat [69], which is predicted from the free energy computed above (recall that quantum integrability and the factorized scattering of the quantum solitons is central to the entire TBA computation). This measurement also provides the mass of the quantum soliton. More details on the specific-heat measurement proposal is given in Appendix C. Additional transport signatures of quantum integrability can be obtained from the boundary dSG model, which is realized by terminating the lattice shown in Fig. 2(e) with two impurity junctions, one each for the upper and lower parts of the array. The quantum transport properties can be calculated along the lines of Refs. [15,49] starting with our TBA. However, inclusion of the boundary in the TBA involves additional technical complications. We intend to report on this in a future publication.

## V. SUMMARY

To summarize, we have provided analytical predictions of the thermodynamic properties of the dSG model and a QEC for simulating it. Faithful analog simulation of QFTs using QECs open possibilities in investigation of their nonperturbative properties. Radio-frequency measurements possible with QECs allow a more detailed verification of the QFTs’ properties. At the same time, QEC simulation of integrable QFTs may be viewed as a method to benchmark these analog simulators by comparing experiments to analytical theoret-

ical predictions. This leads to a systematic investigation of more general QFTs, by systematically including perturbations which break integrability. Using the robust controllability and customizability of QECs, this can be achieved as described below. A massless free-boson QFT can be included by taking the SG circuit after removing the Josephson junctions on the vertical links [upper or lower half of Fig. 2(e)]. A massive free-boson QFT is obtained from the latter by adding linear inductances to the ground at each vertical link. Interactions can be realized by either inductive coupling or capacitive coupling or through Josephson junctions. For instance, interacting QFTs may be obtained by specific arrangements of Josephson junctions in tailored geometries [70–72]. Fermionic modes may potentially also be included using the topological Josephson effect [73]. Finally, the QECs proposed in this work can be generalized to simulate higher dimensional QFTs (integrable and nonintegrable) [74].

## ACKNOWLEDGMENTS

A.R. acknowledges fruitful discussions with J. Belletete, M. Devoret, H. Grabert, M. Pawelkiewicz, and A. D. Stone. A.R. acknowledges the support of the Alexander von Humboldt Foundation. H.S. was supported by the ERC Advanced Grant NuQFT. H.S. thanks F. Lesage and P. Simonetti for (early) collaboration on this project, and T. Albash for discussions.

## APPENDIX A: DOUBLE SINE-GORDON MODEL AT ZERO TEMPERATURE

In this section we consider the dSG model in the presence of external magnetic fields $h_i$’s, which couple to the conserved charges $Q_i$’s of the solitons. First, we calculate the ground-state energy in the limit $h_i \to \infty$ and then consider deviations from this limit using CPT. Subsequently, the results are compared with Bethe ansatz calculations. The charges, $Q_i$’s, of the solitons are given by

$$Q_i = \int_{-\infty}^{\infty} dx \frac{1}{\sqrt{\pi}} \partial_1 \varphi_i, \quad i = 1, 2, \qquad (A1)$$

where $\varphi_{1,2} = (\alpha_1\varphi_1 \pm \alpha_2\varphi_2)/2\sqrt{\pi}$. The coupling to the magnetic field is given by the following action:

$$S_{\text{mag}} = -\frac{h_1}{\sqrt{\pi}} \int d^2x \partial_1 \varphi_1 - \frac{h_2}{\sqrt{\pi}} \int d^2x \partial_1 \varphi_2. \qquad (A2)$$

In limit $h_1/M_0, h_2/M_0 \to \infty$, the interaction term $S_{\text{int}}$ can be neglected. The ground state energy can be computed by completing squares in a Gaussian integral and is given by

$$\mathcal{E} = - \frac{(h_1^2 + h_2^2)A - 2h_1h_2B}{2\pi(A^2 - B^2)}. \qquad (A3)$$

The equilibrium values of the charges $Q_{1,2} = \partial\mathcal{E}/\partial h_{1,2}$ are given by

$$Q_1 = \frac{h_1 A - h_2 B}{\pi (A^2 - B^2)}, \quad Q_2 = \frac{h_2 A - h_1 B}{\pi (A^2 - B^2)}. \qquad (A4)$$

We will be considering the situation when the fields $h_1, h_2$ are so chosen that $Q_2 = 0$ which implies, from the above
equation, $h_2 = Bh_1/A$. Then, the charge $Q_1$ is given by

$$Q_1 = \frac{h_1}{\pi A} = \frac{h_1 \alpha_1^2 \alpha_2^2}{\pi^2 (\alpha_1^2 + \alpha_2^2)} \quad (A5)$$

and the ground-state energy is

$$\mathcal{E}_c = - \frac{h_1^2 \alpha_1^2 \alpha_2^2}{2\pi^2 (\alpha_1^2 + \alpha_2^2)}. \quad (A6)$$

Next, we look at finite $h_1/M_0$, $h_2/M_0$ and consider fluctuations of the fields $\phi_{1,2}$ around their asymptotic values given by Eq. (A4). This is where our results differ from those obtained in Ref. [48]. As will be shown below, the corrections to the limit $h_1/M_0, h_2/M_0 \to \infty$ occur at the second order and lead to logarithmic corrections, while the same for the Fateev model occur at fourth order and give rise to polynomial corrections. We define

$$\phi_i(x) = \sqrt{\pi} Q_i x_1 + \tilde{\phi}_i(x), \quad i = 1, 2, \quad (A7)$$

where $\langle \tilde{\phi}_i \rangle = 0$. The original fields $\phi_{1,2}$ can be written as

$$\phi_{1,2} = \frac{\pi}{\alpha_1} (Q_1 \pm Q_2)x_1 + \tilde{\phi}_{1,2}, \quad (A8)$$

where we have defined $\tilde{\phi}_{1,2} = \sqrt{\pi}(\tilde{\phi}_1 \pm \tilde{\phi}_2)/\alpha_1$. Rewriting the action in terms of $\tilde{\phi}_{1,2}$ and considering the case when $Q_2 = 0$, we have $S_0 + S_{\text{mag}} = \tilde{S}_0 + LR\mathcal{E}_c$, where

$$\tilde{S}_0 = \frac{1}{2} \int d^2x \{(\partial_\mu \tilde{\phi}_1)^2 + (\partial_\mu \tilde{\phi}_2)^2\}, \quad (A9) \\ \mathcal{E}_c = - \frac{h_1^2 \alpha_1^2 \alpha_2^2}{2\pi^2 (\alpha_1^2 + \alpha_2^2)}.$$

Here, $L$ is the spatial extent of the one-dimensional (1D) system and $R$ is the inverse temperature. In the following, we will keep the fluctuating fields $\tilde{\phi}_{1,2}$ and evaluate the corrections due to the interaction term $S_{\text{int}}$ perturbatively. We consider the case when $Q_2 = 0$. The interaction term can be written as

$$S_{\text{int}} = \frac{2M_0}{\pi} \int d^2x \cos(\pi Q_1 x_1 + \alpha_1 \tilde{\phi}_1) \cos(\pi Q_1 x_1 + \alpha_2 \tilde{\phi}_2).$$

This leads to

$$\ln Z = -\mathcal{E}_c LR + \ln Z_0 + \frac{1}{2} \langle S_{\text{int}}^2 \rangle_0, \quad (A10)$$

where we have used the fact that $\langle S_{\text{int}} \rangle_0 = 0$, $Z_0 = \int D\hat{\phi}_1 D\hat{\phi}_2 \exp(-\tilde{S}_0)$ and the averages $\langle \rangle_0$ are with respect to $Z_0$. In order to calculate the correction term, we need the following formula for averages of the vertex operators [75]:

$$\langle e^{i\alpha\phi(x_1)} \dots e^{i\alpha\phi(x_n)} e^{-i\alpha\phi(y_1)} \dots e^{-i\alpha\phi(y_n)} \rangle_0 = \frac{\prod_{i>j}^{n} (|x_i - x_j| |y_1 - y_j|)^{\alpha_1^2/2\pi}}{\prod_{i,j=1}^{n} |x_i - y_j|^{\alpha_2^2/2\pi}}, \quad (A11)$$

where the average is taken with respect to the free-field action. Performing the integral, on the integrable manifold, $\alpha_1^2 + \alpha_2^2 = 4\pi$, we get

$$\langle S_{\text{int}}^2 \rangle_0 = \frac{M_0^2}{2\pi} LR \left\{ -2\gamma_E - 2\ln(2\pi|Q_1|) + 2\ln\left(\frac{R}{2a}\right) \right\},$$

where $\gamma_E$ is the Euler constant and $a$ is a lattice cutoff. Using this result, the specific ground-state energy, up to second order in $M_0$, is given by

$$\mathcal{E} = -\frac{1}{L} \frac{\partial \ln Z}{\partial R} = \mathcal{E}_c + \frac{M_0^2}{2\pi} \ln(2\pi|Q_1|) + \mathcal{E}_0, \quad (A12)$$

where

$$\mathcal{E}_0 = -\frac{1}{L} \frac{\partial \ln Z_0}{\partial R} + \frac{M_0^2}{2\pi} \{\gamma_E + 1 + \ln(R/2a)\}. \quad (A13)$$

Next, we calculate the ground-state energy using Bethe ansatz. Consider the case when the fields are chosen so that $Q_2 = 0$. Due to the presence of the magnetic field, the solitons of the $\phi_1$ fields acquire an additional energy $-h_1Q_1$ and as long as $h_1$ is larger than the mass of the $\phi_1$ solitons, these solitons fill up all possible states within a “Fermi interval” $-B < \theta < B$, where $\theta$ is the rapidity parameter and $B$ is the “Fermi wave vector.” The specific ground-state energy is then given by [59,76,77]

$$\mathcal{E}_c = -\frac{m_s}{2\pi} \int_{-B}^{B} d\theta \cosh(\theta) \epsilon(\theta), \quad (A14)$$

where $m_s$ is the mass of the $\phi_1$ soliton and the quasiparticle energies $\epsilon(\theta)$ satisfy

$$\int_{-B}^{B} d\theta \tilde{K}_c(\theta - \theta') \epsilon(\theta') = h_1 - m_s \cosh(\theta), \quad (A15)$$

with the boundary condition $\epsilon(\pm B) = 0$. Here,

$$\tilde{K}_c(\theta) = \delta(\theta) - \frac{1}{2\pi} \frac{d}{d\theta} [\delta_{p_1}(\theta) + \delta_{p_2}(\theta)], \quad (A16)$$

where $\delta_{p_i}$'s are soliton-soliton scattering phase shifts (for instance, see Ref. [48] for explicit forms). The Fourier transform of $\tilde{K}_c$ is given by

$$K_c(\omega) = \frac{\sinh[\pi\omega(p_1 + p_2)/2]}{2\sinh(\pi\omega p_1/2)\sinh(\pi\omega p_2/2)} \tanh\left(\frac{\pi\omega}{2}\right). \quad (A17)$$

The leading-order energy contribution in the limit of $h_1/M_0, h_2/M_0 \to \infty$ can be obtained by the standard manipulations [48,52] and the result is identical to that of the Fateev model. This leads to

$$\mathcal{E}_c = \frac{h_1^2}{2\pi K_c(0)} = -\frac{h_1^2 p_1 p_2}{\pi(p_1 + p_2)}. \quad (A18)$$

Comparing Eqs. (A6) and (A18), we get

$$p_i = \frac{\alpha_i^2}{2\pi}, \quad i = 1, 2. \quad (A19)$$

In terms of $p_i$'s, the integration manifold is then given by $p_1 + p_2 = 2$. Next, we consider deviations from the $h_1/M_0, h_2/M_0 \to \infty$ limit. This is done by a Wiener-Hopf calculation. Since we are interested in the dSG model on the integration manifold, we set $p_2 = 2 - p_1$ at the outset. Then, the scattering kernel reduces to

$$K_c(\omega) = \frac{\sinh(\pi\omega)}{2\sinh\left(\frac{\pi\omega p_1}{2}\right)\sinh\left(\pi\omega - \frac{\pi\omega p_1}{2}\right)} \tanh\left(\frac{\pi\omega}{2}\right). \quad (A20)$$

The kernel factorizes as $K_c(\omega) = 1/N(\omega)N(-\omega)$ [48], where

$$N(\omega) = \sqrt{\frac{4\pi}{p_1(2-p_1)} \frac{p(i\omega)p\left(\frac{i\omega}{2}\right)e^{i\omega\Lambda}}{p\left(\frac{i\omega p_1}{2}\right)p(i\omega - \frac{i\omega p_1}{2})p\left(\frac{1}{2} + \frac{i\omega}{2}\right)}},$$
where

$$ \Lambda = -\ln 2 + \frac{2 - p_1}{2} \ln(2 - p_1) + \frac{p_1}{2} \ln p_1. \quad (A21) $$

Then, the Bethe ansatz integral equation is reduced to the following linear integral equation for the $v(\omega)$, given by

$$ v(k) = -\frac{ih_1 N(0)}{k} + \frac{im_s e^B N(-i)}{2} + \int_{C_+} \frac{d\omega}{2\pi i} \frac{e^{2i\omega B}}{w+k} \alpha(\omega)v(\omega), \quad (A22) $$

where $m_s$ is the mass of physical excitations of solitons of $\varphi_1$, $\alpha(\omega) = N(\omega)/N(-\omega)$, and

$$ ih_1 N(0) = -\frac{im_s e^B}{2} N(-i) + \int_{C_+} \frac{d\omega}{2\pi i} e^{2i\omega B} \alpha(\omega)v(\omega). \quad (A23) $$

Here, $C_+$ includes all the poles on the positive imaginary axis, but not the one at $\omega = -k$. Then, the ground-state energy is given by

$$ \begin{aligned} \mathcal{E}(h_1) - \mathcal{E}(0) &= -\frac{h_1^2 N(0)^2}{2\pi} - \frac{h_1 N(0)}{\pi} \int_{C_+} \frac{d\omega}{2\pi i} \frac{e^{2i\omega B}}{\omega - i} \alpha(\omega)v(\omega), \end{aligned} $$

where $C'_+$ is the contour that includes the poles included by $C_+$ plus the one at $\omega = i$. Our goal is to calculate this ground-state energy, which gets contribution from the poles of $\alpha(\omega)$, $v(\omega)$, and $\omega = i$. We include these contributions iteratively. Keeping to the lowest order contribution, we evaluate the contribution due to the pole at $\omega = \omega_1 = i$. This leads to

$$ \begin{aligned} \mathcal{E}(h_1) - \mathcal{E}(0) &= -\frac{h_1^2 N(0)^2}{2\pi} + \frac{h_1^2 N(0)^2}{\pi} \int_{C_+} \frac{d\omega}{2\pi i} \frac{e^{2i\omega B}}{\omega(\omega - i)^3} \tilde{\alpha}(\omega), \end{aligned} $$

where $\tilde{\alpha}(\omega) = (\omega - i)\alpha(\omega)$. Application of the residue theorem leads to

$$ \begin{aligned} \mathcal{E}(h_1) - \mathcal{E}(0) &= -\frac{h_1^2 N(0)^2}{2\pi} - \frac{h_1^2 N(0)^2}{2\pi} \frac{m_s^2 \sin^2(\pi p_1/2)}{2h_1^2(2-p_1)p_1} \\ & \times \left\{ -1 + 2\gamma_E + \ln 4 + \frac{4}{p_1-2} - 2\ln \frac{h_1(2-p_1)p_1}{m_s} \right. \\ & \left. + \pi \cot \frac{\pi p_1}{2} + 2\ln \sin \frac{\pi p_1}{2} + 2\psi\left(\frac{p_1}{2}-1\right) \right\}, \quad (A24) \end{aligned} $$

where $\psi$ is the digamma function and $N(0) = \sqrt{(2-p_1)p_1}$. Then, the ground-state energy is given by

$$ \begin{aligned} \mathcal{E}(h_1) - \mathcal{E}(0) &= -\frac{h_1^2 p_1 (2 - p_1)}{2\pi} - \frac{m_s^2 \sin^2(\pi p_1/2)}{4\pi} \\ & \times \left\{ -1 + 2\gamma_E + \ln 4 + \frac{4}{p_1 - 2} \right. \\ & \left. - 2 \ln \frac{h_1(2 - p_1)p_1}{m_s} + \pi \cot \frac{\pi p_1}{2} \right. \\ & \left. + 2 \ln \sin \frac{\pi p_1}{2} + 2 \psi \left( \frac{p_1}{2} - 1 \right) \right\}. \quad (A25) \end{aligned} $$

The contributions to the ground-state energy can be identified as follows. The first term on the right-hand side is $\mathcal{E}_c$, the leading order contribution to the ground-state energy when $h_1/M_0, h_2/M_0 \to \infty$. In the second term, the $h_1$-independent term can be identified as the ground-state energy $\mathcal{E}(0)$, while the logarithmic term dependent on $h_1$ is the perturbative correction as we move away from the infinitely large magnetic fields limit. Thus, we get

$$ \mathcal{E}(h_1) = \mathcal{E}_c + \delta\mathcal{E}, \quad (A26) $$

where

$$ \delta\mathcal{E} = \frac{m_s^2}{2\pi} \sin^2(\pi p_1/2) \ln h_1 \quad (A27) $$

and the ground-state energy in the absence of magnetic field as

$$ \begin{aligned} \mathcal{E}(0) = & \frac{m_s^2 \sin^2(\pi p_1/2)}{4\pi} \left\{ -1 + 2\gamma_E + \ln 4 + \frac{4}{p_1 - 2} \right. \\ & \left. - 2 \ln \frac{(2 - p_1)p_1}{m_s} + \pi \cot \frac{\pi p_1}{2} \right. \\ & \left. + 2 \ln \sin \frac{\pi p_1}{2} + 2\psi\left(\frac{p_1}{2} - 1\right) \right\}. \end{aligned} \quad (A28) $$

Furthermore, comparing the coefficient of $\ln h_1$ term in Eqs. (A12) and (A27) leads to an exact expression connecting the interaction parameter $M_0$ of the action with the masses of the $\varphi_1$ solitons:

$$ M_0 = m_s \sin \frac{\pi p_1}{2}. \quad (A29) $$

Here, we have used the definition of $Q_1$ in Eq. (A5) and the relation Eq. (A19). We note that the masses of the $\varphi_2$ solitons are the same as that of the $\varphi_1$ solitons on the integrability manifold.

## APPENDIX B: DOUBLE SINE-GORDON MODEL AT NONZERO TEMPERATURE

In this section, we analyze the dSG model at finite temperature. The solitons of this model scatter in a factorized manner. In particular, the scattering of the two quantum numbers of each soliton also occurs independent of each other: $S = S_{p_1} \otimes S_{p_2}$, where on the integrable manifold $p_1 + p_2 = 2$. This relation between $p_1, p_2$ forces one of the amplitude to be in the attractive regime and the other to be in the repulsive regime. This is precisely why the analysis is different from the Fateev model, for which the TBA analysis was done only for the case $p_1, p_2 \ge 1$ [48]. Without loss of generality, we choose $p_1 < 1$, i.e., attractive, and $p_2 > 1$ repulsive. In addition to the solitons, the spectrum of the theory also includes $n < 1/p_1$ particles which are neutral with respect to first quantum number, but still carry the second quantum number. They are the bound states of the fundamental “double” solitons and correspond to the ordinary breather poles of the attractive sine-Gordon amplitude of the first sector. The masses of these breathers are given by the standard formula $m_i = 2m_s \sin[\pi i/2(n+1)]$, where $i = 1, \dots, n$ [28]. The TBA analysis of this problem is in general highly nontrivial, scattering in each sector being non-diagonal. We will restrict therefore ourselves to the simpler case $p_1 = 1/(n+1)$ for which the sector one scattering
is diagonal. Our TBA analysis will therefore consist in the diagonalization problem of the second sector only, which can be done by means of the Takahashi-Suzuki classification of the solution of the Bethe equations [56–58].

The periodicity condition for the wave function of this given number of solitons, antisolitons, and breathers on a circle of length $L$ becomes the following set of equations:

$$ 
\begin{aligned}
& \prod_{\delta=s,a} \prod_{k_\delta} S_{p\delta}^{p_1}(\theta - \theta_{k_\delta}) \prod_{i=1,...,n} \prod_{k_{b_i}} S_{pb_i}^{p_1}(\theta - \theta_{k_{b_i}}) \lambda_p^{p_2}(\theta) \\
& = e^{-im_p L \sinh \theta}, \quad p = s, a, \\
& \prod_{\delta=s,a} \prod_{k_\delta} S_{b_i\delta}^{p_1}(\theta - \theta_{k_\delta}) \prod_{j=1,...,n} \prod_{k_{b_j}} S_{b_ib_j}^{p_1}(\theta - \theta_{k_{b_j}}) \lambda_{b_i}^{p_2}(\theta) \\
& = e^{-im_{b_i} L \sinh \theta},
\end{aligned}
\quad (B1) $$

where $i = 1, ..., n$. In this equation, $\theta$ is the rapidity of the particle going around the circle, and $\theta_{k_s}$, $\theta_{k_a}$, $\theta_{k_{b_i}}$ are the rapidities of the incoming soliton, antisoliton, and the breathers respectively. The scattering coefficients for the sine-Gordon model, $S^{p_1}$'s, are well known [28], while $\lambda^{p_2}$ is the contribution of the second phase shift due to the scattering of the second quantum number. This contribution is computed by diagonalization of the transfer matrix of a second sector soliton going around the world repulsively interacting with solitons and antisolitons with amplitudes $S^{p_2}$. The diagonalization procedure for repulsive sine-Gordon model produces in the thermodynamic limit equations for the densities of states in terms of massless pseudoparticles and a physical particle carrying the mass of the soliton. When “gluing” this result with the above equations, the only care needed will be in the identification of the role of the massive physical particles. This is done using the algebraic Bethe ansatz technique [57, 59, 60]. The calculation is nontrivial. First, we do it for $n = 2$ and then generalize the results so obtained. The detailed derivation for $n = 2$ is given at the end of the paper (Appendix D) and here, we present directly the final results:

$$ 
\begin{aligned}
& \epsilon_p(\theta) = \hat{\phi}_n \star \mathcal{L}_{b_n}, \quad p = s, a, \\
& \epsilon_{b_n}(\theta) = \hat{\phi}_n \star (\mathcal{L}_s + \mathcal{L}_a + \mathcal{L}_{b_{n-1}}), \\
& \epsilon_{b_j}(\theta) = \hat{\phi}_n \star (\mathcal{L}_{b_{j+1}} + \mathcal{L}_{b_{j-1}}), \quad j = n-1, \dots, 3, \\
& \epsilon_{b_2}(\theta) = \hat{\phi}_n \star (\mathcal{L}_{b_3} + \mathcal{L}_{b_1} - \mathcal{L}_1), \\
& \epsilon_{b_1}(\theta) = \hat{\phi}_n \star (\mathcal{L}_{b_2} - \mathcal{L}_2), \\
& \epsilon_1(\theta) = -m_{b_1} R \cosh \theta + \hat{\phi}_n \star (\mathcal{L}_{b_2} - \mathcal{L}_2), \\
& \epsilon_2(\theta) = -\hat{\phi}_n \star (\mathcal{L}_{b_1} - \mathcal{L}_1 + \mathcal{L}_3), \\
& \epsilon_j(\theta) = -\hat{\phi}_n \star (\mathcal{L}_{j+1} + \mathcal{L}_{j-1}), \quad j = n-1, \dots, 3, \\
& \epsilon_n(\theta) = -\hat{\phi}_n \star (\mathcal{L}_{n-1} + \mathcal{L}_{n+1} + \mathcal{L}_{n+2}), \\
& \epsilon_k(\theta) = -\hat{\phi}_n \star \mathcal{L}_n, \quad k = n+1, n+2.
\end{aligned}
\quad (B2) $$

In the above equation, the TBA kernel is given by $\hat{\phi}_n(\theta) = (n+1)/\cosh[(n+1)\theta]$, $\mathcal{L}(\theta) = \ln[1+e^{\epsilon(\theta)}]$, and $L(\theta) = \ln[1+e^{-\epsilon(\theta)}]$. Finally, we use the following convention for the convolution: $a \star b = \int \frac{d\theta'}{2\pi} a(\theta - \theta')b(\theta')$. The pictorial representation for the TBA equations is given in Fig. 1. The solid circles denote the physical massive particles, while the empty ones denote the pseudoparticles. The cross on the

empty circle encodes which of the pseudoparticles has a massive term in its TBA equation. The connectivity of the diagram encodes the structure of the TBA equations. Finally, the arrows denote, for the TBA equation for a given particle, the sign of the contribution of the different particles.

Note the somewhat special structure of the TBA equations. In contrast to the usual structure of the TBA equations, in Eq. (B2), only one pseudoparticle (1) has a (negative!) mass term, while the physical massive particles do not. This is because the equations are written in terms of both the functions $\mathcal{L}$ and $L$, which allows us to write the TBA equations concisely in terms of a single universal kernel $\hat{\phi}_n$. When written only in terms of the $L$'s, the TBA equations are in the “standard form” and do have the expected structure.

Note the difference of the TBA equations from those obtained for the Fateev model. The latter model involves, in addition to the two bosonic fields considered in this work, a noncompact bosonic field $\varphi$ with coupling $\beta$. The TBA for the Fateev model has been done in the regime when both the couplings $\alpha_1, \alpha_2$ are larger than $\sqrt{2\pi}$. In this regime, both the sine-Gordon models are in the repulsive regime (without any breathers in either sector) and the resultant TBA is obtained by gluing two sine-Gordon TBA diagrams in the repulsive regime. In contrast, in the dSG model, the pure quantum integrable manifold corresponds to $\alpha_1^2 + \alpha_2^2 = 4\pi$, with $\alpha_1 \neq \alpha_2$. Thus, the breathers in one of the sectors have to be taken into account to correctly obtain the thermodynamics of the dSG model. This difference is visible by taking the naïve limit of $\beta = 0$ in Eq. (33) of Ref. [48], which gives divergent results for the ground-state energy.

To check the TBA equations, we first compute the central charge in the UV limit. The central charge is given by

$$ c_{\text{UV}} = \frac{6}{\pi^2} \int d\theta \cosh \theta \left[ L_s(\theta) + L_a(\theta) + \sum_{j=1}^{n} \frac{m_{b_j}}{m_s} L_{b_j}(\theta) \right]. \quad (B3) $$

For TBA equations with off-diagonal scattering, the central charge can be computed in the usual manner using the solutions of the energies $\epsilon_s, \epsilon_a, \epsilon_{b_i}, \epsilon_i$'s in the UV ($T \to \infty$), and IR ($T \to 0$) limits [61,78]. The resultant expression for the central charge can be written as

$$ c_{\text{UV}} = \frac{6}{\pi^2} \left\{ 2\mathcal{L}_{\text{dlog}}\left(\frac{x_s}{1+x_s}\right) + \sum_{j=2}^{n} \left[ \mathcal{L}_{\text{dlog}}\left(\frac{x_j}{1+x_j}\right) \right. \right. \\
\left. \left. + \mathcal{L}_{\text{dlog}}\left(\frac{x_{b_j}}{1+x_{b_j}}\right) \right] + 2\mathcal{L}_{\text{dlog}}\left(\frac{x_{n+1}}{1+x_{n+1}}\right) \right. \\
\left. + \mathcal{L}_{\text{dlog}}\left(\frac{x_{b_1}}{1+x_{b_1}}\right) + \mathcal{L}_{\text{dlog}}\left(\frac{1}{1+x_1}\right) \right. \\
\left. - \sum_{j=2}^{n} \mathcal{L}_{\text{dlog}}\left(\frac{y_j}{1+y_j}\right) - \mathcal{L}_{\text{dlog}}\left(\frac{1}{1+y_1}\right) \right. \\
\left. - 2\mathcal{L}_{\text{dlog}}\left(\frac{y_{n+1}}{1+y_{n+1}}\right) \right\}, \quad (B4) $$
where $\mathcal{L}_{\text{dlog}}$ is the Rogers dilogarithm function and we have used $\epsilon_s = \epsilon_a, \epsilon_{n+1} = \epsilon_{n+2}$. Furthermore, $x_p = e^{-\epsilon_p}, p = s, a, b_i, i, i = 1, \dots, n$ in the UV limits, while the $y_p = e^{-\epsilon_p}$ in the IR limit. Note that in the IR limit, the energies of the massive physical particles diverge, so they do not contribute to the expression of the central charge [61]. The solutions of the TBA equations in the UV limit are given by

$$ 
\begin{aligned}
x_{n+1} &= \frac{1}{x_s} = n, \quad x_1 = x_{b_1} = 1, \\
x_j &= \frac{1}{x_{b_j}} = j^2 - 1, \quad j = 2, \dots, n,
\end{aligned} 
\quad (B5) 
$$

while for the IR limit, the nonzero y's are given by

$$ 
\begin{aligned}
y_{n+1} &= n - \frac{1}{2}, \quad y_1 = 3, \\
y_j &= \left(j - \frac{1}{2}\right)^2 - 1, \quad j = 2, \dots, n.
\end{aligned} 
\quad (B6) 
$$

Plugging these solutions into Eq. (B4) and using some remarkable dilogarithm identities [79], we arrive at the desired result of

$$ c_{\text{UV}} = 2. \quad (B7) $$

Next, we compute the effective central charge as we move away from the UV limit. The perturbation $S_{\text{int}}$ has dimension $(\alpha_1^2 + \alpha_2^2)/(4\pi)$. Thus, on the integrable manifold, the perturbation has dimension 1. This leads to logarithmic corrections to the central charge [61,62] as one moves away from the UV limit, in addition to the polynomial corrections in powers of $(m_s R)^2$ [54]. We compute these corrections by solving the TBA equations numerically and doing a numerical fit of the resulting expression for the central charge. The coefficient of the logarithmic correction for the different cases is then verified using CPT. These calculations are described below for $n=2$ and $n=3$.

Consider the case $n=2$, i.e., $p_1 = 1/3$, $p_2 = 5/3$. The first sector spectrum consists of soliton $s$, antisoliton $a$, and two breathers $b_1$ and $b_2$ whose masses are $m_{b_1} = m_s$ and $m_{b_2} = \sqrt{3}m_s$. For the numerical solutions of the TBA equations, it is convenient to use the equations in the standard form, given by

$$ 
\begin{aligned}
\epsilon_p(\theta) &= m_p R \cosh \theta + \hat{\phi}_1 \star (L_s + L_a + L_{b_1} - L_1) \\
&\quad + \hat{\phi}_3 \star L_{b_2}, \quad p = s, a, \\
\epsilon_{b_2}(\theta) &= m_{b_2} R \cosh \theta + \hat{\phi}_3 \star (L_s + L_a + L_{b_1} - L_1) \\
&\quad + 2\hat{\phi}_1 \star L_{b_2}, \\
\epsilon_{b_1}(\theta) &= m_{b_1} R \cosh \theta + \hat{\phi}_1 \star (L_s + L_a + L_{b_1} - L_1) \\
&\quad + \hat{\phi}_3 \star L_{b_2} - \hat{\phi}_2 \star L_2, \\
\epsilon_1(\theta) &= \hat{\phi}_1 \star (L_s + L_a + L_{b_1} - L_1) + \hat{\phi}_3 \star L_{b_2} - \hat{\phi}_2 \star L_2, \\
\epsilon_2(\theta) &= -\hat{\phi}_2 \star (L_{b_1} - L_1 + L_2 + L_3), \\
\epsilon_k(\theta) &= -\hat{\phi}_2 \star L_2, \quad k = 3, 4,
\end{aligned} 
\quad (B8) 
$$

where the kernels are given by

$$ 
\begin{aligned}
\hat{\phi}_1(\theta) &= \frac{3}{2 \cosh(3\theta/2)}, \quad \hat{\phi}_2(\theta) = \frac{3}{\cosh(3\theta)}, \\
\hat{\phi}_3(\theta) &= 3\sqrt{2} \frac{\cosh(3\theta/2)}{\cosh(3\theta)}.
\end{aligned} 
\quad (B9) 
$$

<figure><img src="image_6.png" /><figcaption>FIG. 3. Effective central charge for $n=2$. The solid blue line is obtained by solving the TBA equations given by Eq. (B8). The black dots are obtained by performing a fit to Eq. (B10).</figcaption></figure>

We solve Eq. (B8) numerically and plug it in Eq. (B3) to obtain the effective central charge $c_{\text{eff}}$. The result is shown in Fig. 3. The first six correction terms capturing the deviation from the UV limit is obtained by fitting the numerical data to the following expression:

$$ c_{\text{eff}} = c_0 + a_{\text{ln}}^{(2)} \frac{(m_s R)^2}{(2\pi)^2} \ln(m_s R) + \sum_{j=1}^{5} a_j^{(2)} (m_s R)^{2j}, \quad (B10) $$

where the superscript 2 denotes the value of $n$. This leads to

$$ 
\begin{aligned}
c_0 &= 2, \quad a_{\text{ln}}^{(2)} = 6, \quad a_1^{(2)} = -0.39283, \quad a_2^{(2)} = 0.01051, \\
a_3^{(2)} &= -0.01554, \quad a_4^{(2)} = 0.25742, \quad a_5^{(2)} = -1.54212.
\end{aligned} 
$$

Identical calculations can be done for $n=3$. We only provide the final fit to the effective central charge:

$$ 
\begin{aligned}
c_0 &= 2, \quad a_{\text{ln}}^{(3)} = 3.5149, \quad a_1^{(3)} = -0.29775, \quad a_2^{(3)} = 0.00739, \\
a_3^{(3)} &= -0.00818, \quad a_4^{(3)} = 0.1851, \quad a_5^{(3)} = -1.2709.
\end{aligned} 
$$

The coefficient of the logarithmic term in the effective central charge calculation can be checked analytically using CPT. Taking Eq. (A12) for $h_1=0$ and using Eq. (A29), we arrive at the expression for the free energy:

$$ 
\begin{aligned}
f &= -\frac{1}{LR} \ln Z = -\frac{\pi c_{\text{eff}}}{6R^2} \\
&= -\frac{1}{LR} \ln Z_0 - \frac{m_s^2}{\pi} \sin^2 \frac{\pi}{2(n+1)} \ln \left( \frac{R}{2a} \right).
\end{aligned} 
\quad (B11) 
$$

Thus, we arrive at an analytical expression for the logarithmic correction in $c_{\text{eff}}$, proportional to $(m_s R)^2/(2\pi)^2 \ln(m_s R)$, denoted by $a_{\text{ln}}$:

$$ a_{\text{ln}} = 24 \sin^2 \frac{\pi}{2(n+1)}. \quad (B12) $$

For $n=2$ and $n=3$, this coefficient is given by 6 and 3.5147, which are in good agreement with what was obtained by numerically solving the TBA equations [Eqs. (B11) and
<figure><img src="image_7.png" /><figcaption>FIG. 4. Unit cells for a free-boson quantum field theory (a), the quantum sine-Gordon model (b), and the double sine-Gordon model (c).</figcaption></figure>

(B11)]. We have checked our results for n = 4 and n = 5 as well.

## APPENDIX C: DOUBLE SINE-GORDON MODEL WITH QUANTUM CIRCUITS

In this Appendix, we present some additional details on the realization of the dSG with quantum electronic circuits. As mentioned in the main text, the dSG model can be viewed as two sine-Gordon (SG) models coupled by a mutual inductance and a capacitance. Below, we sketch how the different terms of the circuit contribute to the different terms of the action presented in Eq. (1) of the main text.

Consider the circuit in Fig. 4(a). Define the node flux at a point in space-time (x, t) as $\phi(x, t)$ [25], which plays the role of the bosonic field. As shown in Refs. [67,80–82], a 1D array built out of this unit cell, in the continuum limit, provides a circuit realization of a 1D Luttinger liquid. The action of the latter is given by $S_{LL} = \frac{1}{2\pi K} \int d^2x(\partial_\mu\phi)^2$, where we have rescaled the space and time axis to set the plasmon velocity to unity. It is essential that $E_J, E_{C_0} \gg E_{C_J}$ so that the phase-slip rates are exponentially suppressed. This ensures that despite the presence of the nonlinear Josephson element on the horizontal links, the resultant action is Gaussian. Addition of a Josephson junction on vertical links adds the Josephson cosine potential to the action, leading to the sine-Gordon action, where $S_{SG} = S_{LL} - E_{J_0} \int d^2x \cos\phi$, where $E_{J_0}$ is the strength of the Josephson nonlinearity. It is this nonlinearity that gives rise to the cosine nonlinearity of the SG action. In the quantum field theory notation, the SG coupling constant $\beta$ is given by $\beta = \sqrt{4\pi K}$. The free-fermion point of the SG model corresponds to $\beta = \sqrt{4\pi}$, while $\beta \le (\ge) \sqrt{4\pi}$ corresponds to the attractive (repulsive) regime of the model.

From the above circuits, it is straightforward to infer the circuit realizing the dSG action, presented in Fig. 4(c). Consider the upper and lower halves of the array, without the mutual inductance M and the capacitance C, which couples the two halves. From the above explanation, it is clear that each of the two halves of the array realize the SG actions, given by $S_{SG,i}$, $i = 1, 2$:

$$S_{SG,i} = \frac{1}{2\pi K} \int d^2x (\partial_\mu \varphi_i)^2 - E_{J_0} \int d^2x \cos \varphi_i, \quad i = 1, 2. \qquad (C1)$$

Now, consider the mutual inductance M and the capacitance C. The first leads to a term in the action

$$S_M = -M \int d^2x \partial_x \varphi_1 \partial_x \varphi_2, \qquad (C2)$$

while the second leads to

$$S_C = \frac{C}{2} \int d^2x (\partial_t \varphi_1 - \partial_t \varphi_2)^2. \qquad (C3)$$

The last two equations follow directly from the current-voltage constitutive relations of the circuit elements under consideration. The total action of the dSG circuit is given by $S_{SG,1} + S_{SG,2} + S_M + S_C$. From this action, simple algebraic manipulations give rise to the dSG action in terms of the fields $\varphi_{1,2}$, given by

$$S_{dSG} = \frac{A}{2} \int d^2x [(\partial_\mu \varphi_1)^2 + (\partial_\mu \varphi_2)^2] + B \int d^2x \partial_\mu \varphi_1 \partial_\mu \varphi_2 \\ - \frac{M_0}{\pi} \int d^2x \{\cos(2\sqrt{\pi}\varphi_1) + \cos(2\sqrt{\pi}\varphi_2)\}, \qquad (C4)$$

where $A = \pi(1/\alpha_1^2 + 1/\alpha_2^2)$ and $B = \pi(1/\alpha_1^2 - 1/\alpha_2^2)$. In terms of the circuit components,

$$\alpha_1^2 = \frac{K\pi}{2}, \quad \alpha_2^2 = \frac{\pi}{2K + 4\pi C}, \qquad (C5)$$

which are the expressions provided in the main text of the paper.

In the above derivation of the dSG action from the circuit, we have implicitly assumed that the unit cell is repeated an infinite number of times with periodic/open boundary conditions. However, real experiments are done with a finite number of such unit cells. There are two main difficulties in realizing arrays of a very large number of junctions. First, in the derivation of the dSG action, we have neglected the phase slips in the Josephson junctions in the horizontal links. However, for a large enough number of unit cells, this rate is no longer negligible. A crude estimate for the maximum number of unit cells is provided by the inverse of the phase-slip rate and is given by $e^{\sqrt{E_J/E_{C_J}}}$ [64,67]. The same problem also provides a limitation on the frequencies over which the phase-slip rates can be safely neglected and is given by $u\sqrt{C_J/C_0}/a_0$, where $a_0$ is the lattice spacing and $u = a_0\sqrt{2E_{C_0}E_J}$ is the plasmon velocity [67]. Second, the superinductance of the Josephson junction on the horizontal links should not be shunted by the capacitance to ground. This is ensured by choosing $N < C_J/C_0$ [64]. Despite these difficulties, arrays with up to $O(10^5)$ Josephson junctions have been reliably built and experimentally analyzed [83] and we are optimistic that our analytical predictions can be experimentally verified.

Next, we comment on some of the possible measurements in the circuit. Estimation of $\alpha_1, \alpha_2$ can be done by biasing with two current sources (the “magnetic fields” $h_1, h_2$ considered earlier) and measuring the transmission of an input field through the dSG lattice. This provides an estimate of the average current, $\langle\partial_x\varphi_{1,2}\rangle$, flowing through the horizontal links in the circuit. A different approach will be to locally probe the dSG circuit by coupling resonators to a unit cell and measuring current/voltage correlations. While transport properties can be measured in the ways mentioned above,
the specific heat can be measured using AC calorimetry [84]. The latter method involves inducing temperature oscillations on the array patterned on a membrane. The amplitude of the oscillations, together with the applied amplitude and frequency of the applied AC power provides the specific heat [84,85]. A detailed theoretical analysis for the specific-heat measurements for the sine-Gordon model has been made in Refs. [69,86].

We note that opening the circuit to drive and dissipation leads to an interesting and physically relevant problem of dissipative quantum field theories, a largely open field. The calculations made for the sine-Gordon model [87,88] indicate a finite width for the sine-Gordon minima which presumably leads to a finite width for the different particles in the spec-

trum. We speculate that a similar phenomenon may occur in the dSG model. However, a full quantum analysis is still an open problem even for the sine-Gordon model and we hope to return to this problem in the future.

## APPENDIX D: DERIVATION OF THE TBA EQUATIONS

**FOR $n = 2$**

Consider again the case $n = 2$, i.e., $p_1 = 1/3$, $p_2 = 5/3$. The first sector spectrum consists of soliton $s$, antisoliton $a$, and two breathers $b_1$ and $b_2$ whose masses are $m_{b_1} = m_s$ and $m_{b_2} = \sqrt{3}m_s$. We refer to the two sine-Gordon sectors as $T_1$, $T_2$ which correspond to the couplings $\alpha_1, \alpha_2$. The periodicity condition for the wave function of a given number of solitons, antisolitons, and breathers on a circle of length $L$ becomes the following set of equations:

$$ 
\begin{align} 
& \prod_{k_s} S_{ss}^{p_1}(\theta - \theta_{k_s}) \prod_{k_a} S_{sa}^{p_1}(\theta - \theta_{k_a}) \prod_{k_{b_1}} S_{sb_1}^{p_1}(\theta - \theta_{k_{b_1}}) \prod_{k_{b_2}} S_{sb_2}^{p_1}(\theta - \theta_{k_{b_2}}) \lambda_s^{p_2}(\theta | \{\theta_{k_s}\}, \{\theta_{k_a}\}, \{\theta_{k_{b_1}}\}, \{\theta_{k_{b_2}}\}) = e^{-im_s L \sinh \theta}, \nonumber \\
& \prod_{k_s} S_{b_1 s}^{p_1}(\theta - \theta_{k_s}) \prod_{k_a} S_{b_1 a}^{p_1}(\theta - \theta_{k_a}) \prod_{k_{b_1}} S_{b_1 b_1}^{p_1}(\theta - \theta_{k_{b_1}}) \prod_{k_{b_2}} S_{b_1 b_2}^{p_1}(\theta - \theta_{k_{b_2}}) \lambda_{b_1}^{p_2}(\theta | \{\theta_{k_s}\}, \{\theta_{k_a}\}, \{\theta_{k_{b_1}}\}, \{\theta_{k_{b_2}}\}) = e^{-im_1 L \sinh \theta}, \nonumber \\
& \prod_{k_s} S_{b_2 s}^{p_1}(\theta - \theta_{k_s}) \prod_{k_a} S_{b_2 a}^{p_1}(\theta - \theta_{k_a}) \prod_{k_{b_1}} S_{b_2 b_1}^{p_1}(\theta - \theta_{k_{b_1}}) \prod_{k_{b_2}} S_{b_2 b_2}^{p_1}(\theta - \theta_{k_{b_2}}) \lambda_{b_2}^{p_2}(\theta | \{\theta_{k_s}\}, \{\theta_{k_a}\}, \{\theta_{k_{b_1}}\}, \{\theta_{k_{b_2}}\}) = e^{-im_2 L \sinh \theta}, \tag{D1}
\end{align} 
$$

where $\theta$ is the rapidity of the particle going around the world, $\theta_{k_{s,a},b_1,b_2}$ are the rapidities of the external “incoming” solitons, antisolitons, and breathers, and in which $\lambda_i^{p_2}$ represents the contribution to the phase shift of the scattering of the second quantum number. We compute these contributions by the diagonalizing the transfer matrix of a $T_2$ soliton going around the world repulsively interacting with solitons and antisolitons with amplitudes $S_{SG}^{p_2}$. The diagonalization procedure for repulsive sine-Gordon model produces in the thermodynamic limit equations for the densities of states in terms of massless pseudoparticles and a physical particle carrying the mass of the soliton. When gluing this result with the Eq. (D1) the only care needed will be in the identification of the role of the massive physical particles. Through the algebraic Bethe ansatz technique [57,61], the eigenvalues of the repulsive sine-Gordon transfer matrix are

$$ \lambda^{p_2}(\theta, \{\theta_i\}, \{y_r\}) = \prod_i S_{ss}^{p_2}(\theta - \theta_i) \prod_r \frac{\sinh \frac{1}{p_2}(i\pi - y_r - \theta)}{\sinh \frac{1}{p_2}(y_r - \theta)}. \quad (D2) $$

The rapidities $y_r$ label the Bethe eigenvectors of the transfer matrix and have to satisfy the Bethe equation

$$ \prod_i \frac{\sinh \frac{1}{p_2} (i\pi - y_r - \theta)}{\sinh \frac{1}{p_2} (y_r - \theta)} = - \prod_{r'} \frac{\sinh \frac{1}{p_2} (y_r - y_{r'} - i\pi)}{\sinh \frac{1}{p_2} (y_r - y_{r'} + i\pi)}. \quad (D3) $$

The presence of the external particles does not invalidate the string classification [56,58] of the rapidities $y_r$ in the thermodynamic limit. This classification applies actually to the case $p_2 > 2$. It is easy to convince yourself that the classification of the $y_r$'s of Eq. (D3) is equivalent to the classification of the $y_r' = (1 - p_2)y_r$ of an equation like (D3) but with $p_2' = (1 - 1/p_2)^{-1}$. For the general case $p_2 = 2 - 1/(n+1)$ we have the following pattern of $n+2$ strings $y^{(j)}$:

$$ y_k^{(j)} = \text{Re}(y^{(j)}) - \frac{i\pi}{2} \frac{n}{n+1} (2j-1-2k) - \frac{i\pi}{2} p_2 \frac{1-v_j}{2}, \quad k=1,2,\dots,N_j, \quad j=1,\dots,n+2. \quad (D4) $$

The lengths $N_j$ and the parities $v_j$ of the string $y^{(j)}$ are

$$ 
\begin{align} 
N_1 &= 1; \quad N_j = 2j-3, \quad j=2,3,\dots,n+1; \quad N_{n+2} = 2, \nonumber \\ 
v_1 &= 1; \quad v_j = (-1)^{j-1}, \quad j=2,3,\dots,n+1; \quad v_{n+2} = 1. \tag{D5}
\end{align} 
$$

Since there cannot exist different strings with the same length and parity it is useful to denote the strings with $N_v$. The center of a string, i.e., its real part, plays the role in the thermodynamical analysis of the rapidity of a massless pseudoparticle $N_v$ to which we can associate a density $P_{N_v}$ of states containing it and a density of occupied states $P_{N_v}^+$. As usual the Yang-Yang relations between the densities is obtained by plugging the classification (D4) into the Bethe equation (D3) and by taking the imaginary
part of the logarithmic derivative with respect to the external rapidity. The result for the pseudoparticles is

$$P_i(y) = (-1)^{v_i} \left\{ \hat{\phi}_{i0} \star P_0^+ + \sum_j \hat{\phi}_{ij} \star P_j^+ \right\}, \qquad (D6)$$

where

$$\hat{\phi}_{ij}(y) = -id/dy \ln \left\{ \prod_{k_i=1}^{N_i} \prod_{k_j=1}^{N_j} \frac{\sinh \frac{1}{p_2} \left[ y - \frac{i\pi}{2} \frac{n}{n+1} (N_i - N_j - 2k_i + 2k_j) - \frac{i\pi}{2} p_2 \frac{1-v_i v_j}{2} - i\pi \right]}{\sinh \frac{1}{p_2} \left[ y - \frac{i\pi}{2} \frac{n}{n+1} (N_i - N_j - 2k_i + 2k_j) \frac{i\pi}{2} p_2 - \frac{1-v_i v_j}{2} + i\pi \right]} \right\} \\ = (-1)^{v_i v_j} \frac{2}{p_2} \sum_{k_i=1}^{N_i} \sum_{k_j=1}^{N_j} \frac{\sin \frac{\pi}{p_2} \left[ \frac{n}{n+1} (N_i - N_j - 2k_i + 2k_j) + 2 \right]}{\cosh \frac{2y}{p_2} - (-1)^{v_i v_j} \cos \frac{\pi}{p_2} \left[ \frac{n}{n+1} (N_i - N_j - 2k_i + 2k_j) + 2 \right]}, \qquad (D7)$$

$$\hat{\phi}_{i0}(y) = -id/dy \ln \left\{ \prod_{k_i=1}^{N_i} \frac{-\sinh \frac{1}{p_2} \left[ y - \frac{i\pi}{2} \frac{n}{n+1} (N_i + 2 - 2k_i) - \frac{i\pi}{2} p_2 \frac{1-v_i}{2} \right]}{\sinh \frac{1}{p_2} \left[ y - \frac{i\pi}{2} \frac{n}{n+1} (N_i + 2 - 2k_i) \frac{i\pi}{2} p_2 - \frac{1-v_i}{2} - i\pi \right]} \right\} \\ = (-1)^{v_i} \frac{2}{p_2} \sum_{k_i=1}^{N_i} \frac{\sin \frac{\pi}{p_2} \left[ \frac{n}{n+1} (N_i + 2 - 2k_i) \right]}{\cosh \frac{2y}{p_2} - (-1)^{v_i} \cos \frac{\pi}{p_2} \left[ \frac{n}{n+1} (N_i + 2 - 2k_i) \right]} \qquad (D8)$$

are the kernels of the particle i with the physical particle denoted with 0 and with the other pseudoparticles j. The convolution normalization is $a \star b = \int \frac{d\theta}{2\pi} a(\theta - \theta')b(\theta')$. For sine-Gordon model the method would go on with the equation for the density of the physical particle state, obtainable from the periodicity equation. But in the case we are examining now the role of the physical particle is played by the $T_1$ soliton, antisoliton, and breathers. Therefore in our case we have to modify Eq. (D6) to the equation

$$P_i(y) = (-1)^{v_i} \left\{ \phi_{i0} \star (P_s^+ + P_a^+) + \sum_j \hat{\phi}_{ib_j} \star P_{b_j}^+ + \sum_j \hat{\phi}_{ij} \star P_j^+ \right\}, \qquad (D9)$$

where the breather kernels are

$$\hat{\phi}_{ib_j}(\theta) = \hat{\phi}_{i0}\left(\theta + \frac{i\pi}{2} \frac{n-j+1}{n+1}\right) + \hat{\phi}_{i0}\left(\theta - \frac{i\pi}{2} \frac{n-j+1}{n+1}\right). \qquad (D10)$$

The equations for the densities of the physical particles states can be obtained by taking the logarithmic derivative of Eq. (D1) with the explicit formulas for the eigenvalues (D2) and (D4). The result is

$$P_s(\theta) = \frac{m_s L}{2\pi} \cosh \theta + \hat{\phi}_{ss} \star (P_s^+ + P_a^+) + \sum_j \hat{\phi}_{sb_j} \star P_{b_j}^+ + \sum_j \hat{\phi}_{sj} \star P_j^+, \qquad (D11)$$

$$P_{b_i}(\theta) = \frac{m_i L}{2\pi} \cosh \theta + \hat{\phi}_{b_i s} \star (P_s^+ + P_a^+) + \sum_j \hat{\phi}_{b_i j} \star P_{b_i}^+ + \sum_j \hat{\phi}_{b_i j} \star P_j^+. \qquad (D12)$$

The new kernels are

$$\hat{\phi}_{ss}(\theta) = -i \frac{d}{d\theta} \ln S_{ss}^{p_1}(\theta) + \phi_{00}(\theta),$$

$$\hat{\phi}_{sb_j}(\theta) = -i \frac{d}{d\theta} \ln S_{sb_j}^{p_1}(\theta) + \left[ \phi_{00}\left(\theta + \frac{i\pi}{2} \frac{n-j+1}{n+1}\right) + \hat{\phi}_{00}\left(\theta - \frac{i\pi}{2} \frac{n-j+1}{n+1}\right) \right],$$

$$\hat{\phi}_{si}(\theta) = \hat{\phi}_{i0}(\theta),$$

$$\hat{\phi}_{b_i s}(\theta) = \hat{\phi}_{sb_i}(\theta),$$

$$\hat{\phi}_{b_i b_j}(\theta) = -i \frac{d}{d\theta} \ln S_{b_i b_j}^{p_1}(\theta) + \left[ \hat{\phi}_{00}\left(\theta + \frac{i\pi}{2} \frac{2n-i-j+2}{n+1}\right) + \hat{\phi}_{00}\left(\theta + \frac{i\pi}{2} \frac{i-j}{n+1}\right) \right. \\ \left. + \hat{\phi}_{00}\left(\theta - \frac{i\pi}{2} \frac{i-j}{n+1}\right) + \hat{\phi}_{00}\left(\theta - \frac{i\pi}{2} \frac{2n-i-j+2}{n+1}\right) \right],$$

$$\hat{\phi}_{b_i j}(\theta) = \hat{\phi}_{j b_i}(\theta), \qquad (D12)$$
where we denote with

$$ \hat{\phi}_{00}(\theta) = -i \frac{d}{d\theta} \ln S_{ss}^{p_2}(\theta) = \int dk e^{ik\theta} \frac{\sinh(p_2 - 1)^{\frac{\pi k}{2}}}{2 \sinh p_2^{\frac{\pi k}{2}} \cosh \frac{\pi k}{2}} \quad (D13) $$

the soliton-soliton kernel of the pure $p_2$-sine-Gordon TBA. In order to be concrete, let's come back to the specific example $p_1 = 1/3$ and $p_2 = 5/3$. The equations for the density can be written in Fourier space with $\hat{\phi}(\theta) = \int dk e^{ik\theta} \tilde{\phi}(k)$ and assume the form

$$ \begin{aligned} \mathcal{F}T\left[P_s - \frac{m_s L}{2\pi} \cosh \theta\right] &= -\frac{2c_{1/6} s_{1/3}}{s_{5/6}} (\tilde{P}_s^+ + \tilde{P}_a^+) - \frac{s_{1/2}}{s_{5/6}} \tilde{P}_{b_1}^+ - \frac{4c_{1/6}^2 s_{1/3}}{s_{5/6}} \tilde{P}_{b_2}^+ + \frac{s_{1/2}}{s_{5/6}} \tilde{P}_1^+ - \frac{s_{1/3}}{s_{5/6}} \tilde{P}_{1-}^+ + \frac{s_{1/6}}{s_{5/6}} \tilde{P}_2^+ - \frac{s_{1/6}}{s_{5/6}} \tilde{P}_3^+, \\ \mathcal{F}T\left[P_{b_1} - \frac{m_{b_1} L}{2\pi} \cosh \theta\right] &= -\frac{s_{1/2}}{s_{5/6}} (\tilde{P}_s^+ + \tilde{P}_a^+) - \frac{s_{1/6}}{s_{5/6}} \tilde{P}_{b_1}^+ - \frac{2c_{1/6} s_{1/2}}{s_{5/6}} \tilde{P}_{b_2}^+ + \frac{s_{1/6}}{s_{5/6}} \tilde{P}_1^+ - \frac{s_{2/3}}{s_{5/6}} \tilde{P}_{1-}^+ + \frac{2c_{1/3} s_{1/6}}{s_{5/6}} \tilde{P}_2^+ - \frac{2c_{1/3} s_{1/6}}{s_{5/6}} \tilde{P}_3^+, \\ \mathcal{F}T\left[P_{b_2} - \frac{m_{b_2} L}{2\pi} \cosh \theta\right] &= -\frac{4c_{1/6}^2 s_{1/3}}{s_{5/6}} (\tilde{P}_s^+ + \tilde{P}_a^+) - \frac{2c_{1/6} s_{1/2}}{s_{5/6}} \tilde{P}_{b_1}^+ + \frac{c_{1/6}(s_{1/6} + s_{1/3} - 3s_{5/6})}{c_{1/2} s_{5/6}} \tilde{P}_{b_2}^+ + \frac{2c_{1/6} s_{1/2}}{s_{5/6}} \tilde{P}_1^+ \\ &\quad - \frac{2c_{1/6} s_{1/3}}{s_{5/6}} \tilde{P}_{1-}^+ + \frac{s_{1/3}}{s_{5/6}} \tilde{P}_2^+ - \frac{s_{1/3}}{s_{5/6}} \tilde{P}_3^+, \\ &\quad \tilde{P}_1 = \frac{s_{1/2}}{s_{5/6}} (\tilde{P}_s^+ + \tilde{P}_a^+) + \frac{s_{1/6}}{s_{5/6}} \tilde{P}_{b_1}^+ + \frac{2c_{1/6} s_{1/2}}{s_{5/6}} \tilde{P}_{b_2}^+ - \frac{s_{1/6}}{s_{5/6}} \tilde{P}_1^+ + \frac{s_{2/3}}{s_{5/6}} \tilde{P}_{1-}^+ - \frac{2c_{1/3} s_{1/6}}{s_{5/6}} \tilde{P}_2^+ + \frac{2c_{1/3} s_{1/6}}{s_{5/6}} \tilde{P}_3^+, \\ &\quad \tilde{P}_{1-} = \frac{s_{1/3}}{s_{5/6}} (\tilde{P}_s^+ + \tilde{P}_a^+) + \frac{2c_{1/3} s_{1/6}}{s_{5/6}} \tilde{P}_{b_1}^+ + \frac{2c_{1/6} s_{1/3}}{s_{5/6}} \tilde{P}_{b_2}^+ - \frac{s_{2/3}}{s_{5/6}} \tilde{P}_1^+ + \frac{s_{1/6}}{s_{5/6}} \tilde{P}_{1-}^+ + \frac{2c_{1/2} s_{1/6}}{s_{5/6}} \tilde{P}_2^+ - \frac{2c_{1/2} s_{1/6}}{s_{5/6}} \tilde{P}_3^+, \\ &\quad \tilde{P}_2 = \frac{s_{1/6}}{s_{5/6}} (\tilde{P}_s^+ + \tilde{P}_a^+) + \frac{2c_{1/3} s_{1/6}}{s_{5/6}} \tilde{P}_{b_1}^+ + \frac{s_{1/3}}{s_{5/6}} \tilde{P}_{b_2}^+ - \frac{2c_{1/3} s_{1/6}}{s_{5/6}} \tilde{P}_1^+ - \frac{2c_{1/2} s_{1/6}}{s_{5/6}} \tilde{P}_{1-}^+ + \frac{c_{1/2} s_{1/6}}{c_{1/6} s_{5/6}} \tilde{P}_2^+ - \frac{c_{1/2} s_{1/6}}{c_{1/6} s_{5/6}} \tilde{P}_3^+, \\ &\quad \tilde{P}_3 = \tilde{P}_2, \end{aligned} \quad (D14) $$

where $s_m = \sinh m\pi k$ and $c_m = \cosh m\pi k$ and the positive parity of the strings $1_+, 2_+$, and $3_+$ is understood in the notation. Furthermore we never write the corresponding equation for the antisoliton density being trivially $P_s = P_a$.

Simple manipulations with trigonometric identities result in a drastic simplification of the system:

$$ \begin{aligned} \mathcal{F}T\left[P_s - \frac{m_s L}{2\pi} \cosh \theta\right] &= -\frac{1}{2c_{1/3}} (\tilde{P}_s^+ + \tilde{P}_a^+) - \frac{1}{2c_{1/3}} + \tilde{P}_{b_1}^+ - \frac{c_{1/6}}{c_{1/3}} + \tilde{P}_{b_2}^+ - \frac{1}{2c_{1/3}} + \tilde{P}_1^-, \\ \mathcal{F}T\left[P_{b_2} - \frac{m_{b_2} L}{2\pi} \cosh \theta\right] &= -\frac{c_{1/6}}{c_{1/3}} (\tilde{P}_s^+ + \tilde{P}_a^+) - \frac{c_{1/6}}{c_{1/3}} + \tilde{P}_{b_1}^+ - \frac{1}{c_{1/3}} + \tilde{P}_{b_2}^+ - \frac{c_{1/6}}{c_{1/3}} + \tilde{P}_1^-, \\ \mathcal{F}T\left[P_{b_1} - \frac{m_{b_1} L}{2\pi} \cosh \theta\right] &= -\frac{1}{2c_{1/3}} (\tilde{P}_s^+ + \tilde{P}_a^+) - \frac{1}{2c_{1/3}} + \tilde{P}_{b_1}^+ - \frac{c_{1/6}}{c_{1/3}} + \tilde{P}_{b_2}^+ - \frac{1}{2c_{1/3}} + \tilde{P}_1^-, \\ \tilde{P}_1 &= \frac{1}{2c_{1/3}} (\tilde{P}_s^+ + \tilde{P}_a^+) + \frac{1}{2c_{1/3}} + \tilde{P}_{b_1}^+ + \frac{c_{1/6}}{c_{1/3}} + \tilde{P}_{b_2}^+ + \frac{1}{2c_{1/3}} + \tilde{P}_1^- - \frac{1}{2c_{1/6}} + \tilde{P}_{1-}^-, \\ \tilde{P}_{1-} &= \frac{1}{2c_{1/6}} [\tilde{P}_{b_1}^+ + \tilde{P}_1^- + \tilde{P}_2^+ + \tilde{P}_3^+], \\ \tilde{P}_2 &= \frac{1}{2c_{1/6}} \tilde{P}_{1-}^-, \\ \tilde{P}_3 &= \tilde{P}_2. \end{aligned} \quad (D15) $$

The hole densities $P_i^-$ are defined as $P_i^- = P_i - P_i^+$.

We now have all the ingredients to make a thermodynamical analysis. The internal energy is

$$ U = \int d\theta \{m_s \cosh \theta [P_s^+(\theta) + P_a^+(\theta)] + m_{b_1} \cosh \theta P_{b_1}^+(\theta) + m_{b_2} \cosh \theta P_{b_2}^+(\theta)\} \quad (D16) $$

and the entropy, in the Stirling approximation, is

$$ S = \int d\theta \sum_i (P_i \ln P_i + P_i^+ \ln P_i^+ + P_i^- \ln P_i^-). \quad (D17) $$

Next, we minimizing the free energy $F = U - TS$ (here T is the temperature) with respect to the densities given by

$$ \rho_i^+ = \{P_s^+, P_a^+, P_{b_1}^+, P_{b_2}^+, P_3^+, P_1^-, P_{1-}^-\}. \quad (D18) $$
Taking into account the constraints of Eq. (D15), we end up with the following TBA system for the pseudoenergies $\epsilon_i(\theta)$ defined as $\rho_i^+/P_i = e^{-\epsilon_i}/(1 + e^{-\epsilon_i})$:

$$
\begin{align*}
\epsilon_s(\theta) &= m_s R \cosh \theta + \hat{\phi}_1 \star [L_s + L_a + L_{b_1} - L_1] + \hat{\phi}_3 \star L_{b_2}, \\
\epsilon_{b_2}(\theta) &= m_{b_2} R \cosh \theta + \hat{\phi}_3 \star [L_s + L_a + L_{b_1} - L_1] + 2\hat{\phi}_1 \star L_{b_2}, \\
\epsilon_{b_1}(\theta) &= m_{b_1} R \cosh \theta + \hat{\phi}_1 \star [L_s + L_a + L_{b_1} - L_1] + \hat{\phi}_3 \star L_{b_2} - \phi_2 \star L_{1-}, \\
\epsilon_1(\theta) &= \hat{\phi}_1 \star [L_s + L_a + L_{b_1} - L_1] + \hat{\phi}_3 \star L_{b_2} - \hat{\phi}_2 \star L_{1-}, \\
\epsilon_{1-}(\theta) &= -\hat{\phi}_2 \star [L_{b_1} - L_1 + L_2 + L_3], \\
\epsilon_2(\theta) &= -\hat{\phi}_2 \star L_{1-}, \\
\epsilon_3(\theta) &= -\hat{\phi}_2 \star L_{1-},
\end{align*}
\tag{D19}
$$

where R is the inverse temperature. Here we have introduced the notation $L_i(\theta) = \ln(1 + e^{-\epsilon_i})$ and the new kernels $\tilde{\phi}_1(k) = 1/2 \cosh \frac{\pi k}{3}, \tilde{\phi}_2(k) = 1/2 \cosh \frac{\pi k}{6}$, and $\tilde{\phi}_3(k) = \cosh \frac{\pi k}{6} / \cosh \frac{\pi k}{3}$.

A further trivial analytic and algebraic manipulation gives the final form for the $\epsilon$ system:

$$
\begin{align*}
\epsilon_s(\theta) &= \hat{\phi}_2 \star \mathcal{L}_{b_2}, \\
\epsilon_{b_2}(\theta) &= \hat{\phi}_2 \star [\mathcal{L}_s + \mathcal{L}_a + \mathcal{L}_{b_1} - \mathcal{L}_1], \\
\epsilon_{b_1}(\theta) &= \hat{\phi}_2 \star [\mathcal{L}_{b_2} - \mathcal{L}_2], \\
\epsilon_1(\theta) &= -m_{b_1} R \cosh \theta + \hat{\phi}_2 \star [\mathcal{L}_{b_2} - \mathcal{L}_{1-}], \\
\epsilon_{1-}(\theta) &= -\hat{\phi}_2 \star [\mathcal{L}_{b_1} - \mathcal{L}_1 + \mathcal{L}_2 + \mathcal{L}_3], \\
\epsilon_2(\theta) &= -\hat{\phi}_2 \star \mathcal{L}_{1-}, \\
\epsilon_3(\theta) &= -\hat{\phi}_2 \star \mathcal{L}_{1-}.
\end{align*}
\tag{D20}
$$

The functions $\mathcal{L}_i$ are the positive part of $\epsilon_i$, i.e., $\mathcal{L}_i = \ln(1 + e^{\epsilon_i})$. The universal kernel is $\tilde{\phi}_2(k) = 1/2 \cosh \frac{\pi k}{6}$. The above procedure can be generalized to other integer values of $n$, leading to the equations presented in Eq. (B2).

[1] R. P. Feynman, Int. J. Theor. Phys. 21, 467 (1982).

[2] S. Lloyd, Science 273, 1073 (1996).

[3] D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 79, 2586 (1997).

[4] B. Douçot, L. B. Ioffe, and J. Vidal, Phys. Rev. B 69, 214501 (2004).

[5] H. P. Büchler, M. Hermele, S. D. Huber, M. P. A. Fisher, and P. Zoller, Phys. Rev. Lett. 95, 040402 (2005).

[6] J. I. Cirac, P. Maraner, and J. K. Pachos, Phys. Rev. Lett. 105, 190403 (2010).

[7] J. Casanova, L. Lamata, I. L. Egusquiza, R. Gerritsma, C. F. Roos, J. J. García-Ripoll, and E. Solano, Phys. Rev. Lett. 107, 260501 (2011).

[8] S. P. Jordan, K. S. M. Lee, and J. Preskill, Quantum Inf. Comput. 14, 1014 (2014).

[9] S. P. Jordan, K. S. M. Lee, and J. Preskill, Science 336, 1130 (2012).

[10] A. A. Houck, H. E. Türeci, and J. Koch, Nat. Phys. 8, 292 (2012).

[11] M. H. Devoret and R. J. Schoelkopf, Science 339, 1169 (2013).

[12] P. Fendley, A. W. W. Ludwig, and H. Saleur, Phys. Rev. B 52, 8934 (1995).

[13] P. Fendley, A. W. W. Ludwig, and H. Saleur, Phys. Rev. Lett. 74, 3005 (1995).

[14] F. Milliken, C. Umbach, and R. Webb, Solid State Commun. 97, 309 (1996).

[15] F. Lesage, H. Saleur, and P. Simonetti, Phys. Rev. B 56, 7598 (1997).

[16] D. Giuliano and P. Sodano, Nucl. Phys. B 711, 480 (2005).

[17] V. Gritsev, A. Polkovnikov, and E. Demler, Phys. Rev. B 75, 174511 (2007).

[18] A. P. Bukhvostov and L. N. Lipatov, Nucl. Phys. B 180, 116 (1981) [Pisma Zh. Eksp. Teor. Fiz. 31, 138 (1980)].

[19] H. Saleur, J. Phys. A: Math. Gen. 32, L207 (1999).

[20] T. Giamarchi, *Quantum Physics in One Dimension*, International Series of Monographs on Physics (Clarendon, Oxford, 2003).

[21] A. Gogolin, A. Nersesyan, and A. Tsvelik, *Bosonization and Strongly Correlated Systems* (Cambridge University Press, Cambridge, UK, 2004).

[22] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).

[23] S. Sachdev, *Quantum Phase Transitions* (Cambridge University Press, Cambridge, UK, 2011).

[24] M. Tinkham, *Introduction to Superconductivity*, 2nd ed., Dover Books on Physics (Dover, New York, 2004).

[25] M. H. Devoret, *Quantum Fluctuations in Electrical Circuits (Les Houches Session LXIII)* (Elsevier, New York, 1997), pp. 351–386.

[26] This is true also if one tries to “backtrack” to the fermionic formulation, where “non-Fermi liquid” fixed points are then encountered.
[27] D. Iagolnitzer, Phys. Rev. D **18**, 1275 (1978).

[28] A. B. Zamolodchikov and A. B. Zamolodchikov, Ann. Phys. **120**, 253 (1979).

[29] A. Zamolodchikov, Nucl. Phys. B **342**, 695 (1990).

[30] F. Smirnov, *Form Factors in Completely Integrable Models of Quantum Field Theory*, Advanced Series in Mathematical Physics (World Scientific, Singapore, 1992).

[31] T. Senthil and M. P. A. Fisher, Phys. Rev. B **62**, 7850 (2000).

[32] C. Xu and L. Fu, Phys. Rev. B **81**, 134435 (2010).

[33] A. Roy, B. M. Terhal, and F. Hassler, Phys. Rev. Lett. **119**, 180508 (2017).

[34] Y. V. Nazarov and Y. M. Blanter, *Quantum Transport: Introduction to Nanoscience* (Cambridge University Press, Cambridge, UK, 2009).

[35] R. Rajaraman, *Solitons and Instantons: An Introduction to Solitons and Instantons in Quantum Field Theory*, North-Holland personal library (North-Holland, Amsterdam, 1982).

[36] M. Ablowitz and H. Segur, *Solitons and the Inverse Scattering Transform*, SIAM Studies in Applied Mathematics (Society for Industrial and Applied Mathematics, Philadelphia, 2006).

[37] S. Mandelstam, Phys. Rev. D **11**, 3026 (1975).

[38] S. Ghoshal and A. Zamolodchikov, Int. J. Mod. Phys. A **09**, 3841 (1994).

[39] K. Zarembo, Les Houches Summer School: Integrability: From Statistical Systems to Gauge Theory Les Houches, France, June 6-July 1, 2016, Report No. NORDITA 2017-137, UUITP-04/17 (2017), arXiv:1712.07725.

[40] D. Bernard and A. LeClair, Commun. Math. Phys. **142**, 99 (1991).

[41] A. Belavin, A. Polyakov, and A. Zamolodchikov, Nucl. Phys. B **241**, 333 (1984).

[42] A. B. Zamolodchikov, JETP Lett. **46**, 160 (1987) [Pisma Zh. Eksp. Teor. Fiz. **46**, 129 (1987)].

[43] A. Zamolodchikov, in *Integrable Systems of Quantum Field Theory*, edited by M. Jimbo, T. Miwa, and A. Tsuchiya (Academic, San Diego, 1989), pp. 641–674.

[44] P. Francesco, P. Di Francesco, P. Mathieu, D. Sénéchal, and D. Senechal, *Conformal Field Theory*, Graduate Texts in Contemporary Physics (Springer, New York, 1997).

[45] G. Mussardo, *Statistical Field Theory: An Introduction to Exactly Solved Models in Statistical Physics*, Oxford Graduate Texts (Oxford University Press, Oxford, 2010).

[46] I. Affleck, M. Oshikawa, and H. Saleur, Nucl. Phys. B **594**, 535 (2001).

[47] M. Ameduri and C. J. Efthimiou, *J. Nonlinear Math. Phys.* **5**, 132 (1998).

[48] V. A. Fateev, Nucl. Phys. B **473**, 509 (1996).

[49] F. Lesage, H. Saleur, and P. Simonetti, Phys. Rev. B **57**, 4694 (1998).

[50] M. Ameduri, C. J. Efthimiou, and B. Gerganov, Mod. Phys. Lett. A **14**, 2341 (1999).

[51] V. V. Bazhanov, S. L. Lukyanov, and B. A. Runov, Nucl. Phys. B **927**, 468 (2018).

[52] V. Fateev, E. Onofri, and A. Zamolodchikov, Nucl. Phys. B **406**, 521 (1993).

[53] On the integrable manifold, the mass of the $\phi_2$ soliton is also given by $m_s$.

[54] T. R. Klassen and E. Melzer, Nucl. Phys. B **350**, 635 (1991).

[55] R. Dashen, S.-k. Ma, and H. J. Bernstein, Phys. Rev. **187**, 345 (1969).

[56] M. Takahashi and M. Suzuki, Prog. Theor. Phys. **48**, 2187 (1972).

[57] P. Fendley and K. Intriligator, Nucl. Phys. B **372**, 533 (1992).

[58] M. Takahashi, *Thermodynamics of One-Dimensional Solvable Models* (Cambridge University Press, Cambridge, UK, 2005).

[59] V. Korepin, N. Bogoliubov, and A. Izergin, *Quantum Inverse Scattering Method and Correlation Functions*, Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, UK, 1997).

[60] L. D. Faddeev, in *Relativistic Gravitation and Gravitational Radiation, Proceedings, School of Physics, Les Houches, France, September 26-October 6, 1995* (Cambridge University Press, Cambridge, UK, 1997), pp. 149–219.

[61] A. Zamolodchikov, Nucl. Phys. B **358**, 497 (1991).

[62] P. Fendley, H. Saleur, and A. B. Zamolodchikov, Int. J. Mod. Phys. A **08**, 5717 (1993).

[63] A. B. Zamolodchikov, J. Phys. A: Math. Gen. **39**, 12847 (2006).

[64] V. E. Manucharyan, J. Koch, L. I. Glazman, and M. H. Devoret, Science **326**, 113 (2009).

[65] J. Koch, V. Manucharyan, M. H. Devoret, and L. I. Glazman, Phys. Rev. Lett. **103**, 217004 (2009).

[66] K. Le Hur, Phys. Rev. B **85**, 140506(R) (2012).

[67] M. Goldstein, M. H. Devoret, M. Houzet, and L. I. Glazman, Phys. Rev. Lett. **110**, 017002 (2013).

[68] J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A **76**, 042319 (2007).

[69] C. Guarcello, P. Solinas, A. Braggio, and F. Giazotto, Sci. Rep. **8**, 12287 (2018).

[70] N. Bergeal, R. Vijay, V. E. Manucharyan, I. Siddiqi, R. J. Schoelkopf, S. M. Girvin, and M. H. Devoret, Nat. Phys. **6**, 296 (2010).

[71] A. Roy, L. Jiang, A. D. Stone, and M. Devoret, Phys. Rev. Lett. **115**, 150503 (2015).

[72] A four-field integrable QFT with QECs will be described elsewhere.

[73] L. Fu, Phys. Rev. Lett. **104**, 056402 (2010).

[74] A. B. Zamolodchikov, Commun. Math. Phys. **79**, 489 (1981).

[75] A. B. Zamolodchikov, Int. J. Mod. Phys. A **10**, 1125 (1995).

[76] A. Polyakov and P. Wiegmann, Phys. Lett. B **131**, 121 (1983).

[77] P. B. Wiegmann, Phys. Lett. B **152**, 209 (1985).

[78] A. Zamolodchikov, Nucl. Phys. B **358**, 524 (1991).

[79] L. Lewin, *Polylogarithms and Associated Functions* (North-Holland, Amsterdam, 1981).

[80] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Rev. Mod. Phys. **82**, 1155 (2010).

[81] A. Roy and M. Devoret, C. R. Phys. **17**, 740 (2016).

[82] A. Roy and M. Devoret, Phys. Rev. B **98**, 045405 (2018).

[83] R. Kuzmin, N. Grabon, N. Mehta, R. Mencia, N. Pankratova, M. Goldstein, and V. Manucharyan, *Bulletin of American Physical Society* **H33.003** (2018).

[84] F. Fominaya, T. Fournier, P. Gandit, and J. Chaussy, Rev. Sci. Instrum. **68**, 4191 (1997).

[85] H. Rabani, F. Taddei, O. Bourgeois, R. Fazio, and F. Giazotto, Phys. Rev. B **78**, 012503 (2008).

[86] C. Guarcello, P. Solinas, A. Braggio, and F. Giazotto, Phys. Rev. B **98**, 104501 (2018).

[87] K. Hida and U. Eckern, Phys. Rev. B **30**, 4096 (1984).

[88] K. Hida, Z. Phys. B: Condens. Matter **61**, 223 (1985).