Quantum Electronic Circuit Simulation of Generalized sine-Gordon Models

TL;DR

This paper proposes quantum simulators using superconducting circuits to study complex nonlinear quantum field theories, specifically the double sine-Gordon model, and analyzes its thermodynamic properties and quantum integrability.

Contribution

It introduces a novel approach to simulate generalized sine-Gordon models with superconducting circuits and provides analytical and theoretical analysis of their thermodynamic behavior.

Findings

Analytical computation of mass-spectrum and ground state energy.

Derivation of thermodynamic Bethe ansatz equations.

Proposal of experiments to verify quantum integrability.

Abstract

Investigation of strongly interacting, nonlinear quantum field theories (QFT-s) remains one of the outstanding challenges of modern physics. Here, we describe analog quantum simulators for nonlinear QFT-s 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 semi-classical description - a property that is generic to many multi-field QFT-s. 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.

Equations

$${\mathcal{S}}_{\mathrm{dSG}}=\int d^{2}x\frac{1}{2}\sum _{i=1,2}({\partial }_{\mu }{\varphi }_i{)}^2+\frac{2M_0}{\pi }\prod _{i=1,2}\cos ({\alpha }_i{\varphi }_i),$$

$${\cal E}(h_{1})=-\frac{h_{1}^{2}p_{1}(2-p_{1})}{2\pi}+\frac{m_{s}^{2}}{2\pi}\sin^{2}\Big{(}\frac{\pi p_{1}}{2}\Big{)}\ln h_{1}.$$

$${ϵ}_{\gamma }(\theta )$$

$${ϵ}_{b_n}(\theta )$$

$${ϵ}_{b_j}(\theta )$$

$${ϵ}_{b_2}(\theta )$$

$${ϵ}_{b_1}(\theta )$$

$${ϵ}_1(\theta )$$

$${ϵ}_2(\theta )$$

$${ϵ}_j(\theta )$$

$${ϵ}_n(\theta )$$

$${ϵ}_k(\theta )$$

$$(\Delta c_{\rm{eff}})_{\rm{log}}=\frac{6(m_{s}R)^{2}}{\pi^{2}}\sin^{2}\Big{[}\frac{\pi}{2(n+1)}\Big{]}\ln(m_{s}R).$$

$$Q_i={\int }_{-\infty }^{\infty }dx\frac{1}{\sqrt{\pi }}{\partial }_1{\phi }_i,i=1,2,$$

$${\mathcal{S}}_{\mathrm{mag}}$$

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

$$Q_1=\frac{h_{1}A-h_{2}B}{\pi (A^2-B^2)},Q_2=\frac{h_{2}A-h_{1}B}{\pi (A^2-B^2)}.$$

$$Q_1=\frac{h_1}{\pi A}=\frac{h_1{\alpha }_1^2{\alpha }_2^2}{{\pi }^2({\alpha }_1^2+{\alpha }_2^2)}$$

$${\mathcal{E}}_c=-\frac{h_1^2{\alpha }_1^2{\alpha }_2^2}{2{\pi }^2({\alpha }_1^2+{\alpha }_2^2)}.$$

$${\phi }_i(x)=\sqrt{\pi }{Q}_i{x}_1+{\tilde{\phi }}_i(x),i=1,2,$$

$${\varphi }_{1,2}=\frac{\pi }{{\alpha }_1}(Q_1\pm Q_2)x_1+{\tilde{\varphi }}_{1,2},$$

$${\tilde{\mathcal{S}}}_0$$

$${\mathcal{E}}_c$$

$${\mathcal{S}}_{\mathrm{int}}=\frac{2M_0}{\pi }\int d^{2}x\cos (\pi Q_1{x}_1+{\alpha }_1{\tilde{\varphi }}_1)\cos (\pi Q_1{x}_1+{\alpha }_2{\tilde{\varphi }}_2).$$

$$\ln Z$$

$$⟨e^{i\alpha \varphi (x_1)}\dots e^{i\alpha \varphi (x_n)}{e}^{-i\alpha \varphi (y_1)}\dots e^{-i\alpha \varphi (y_n)}{⟩}_0$$

$$\langle{\cal S}_{\rm{int}}^{2}\rangle_{0}=\frac{M_{0}^{2}}{2\pi}LR\Big{\{}-2\gamma_{E}-2\ln(2\pi|Q_{1}|)+2\ln\Big{(}\frac{R}{2a}\Big{)}\Big{\}},$$

$$\mathcal{E}$$

$${\cal E}_{0}=-\frac{1}{L}\frac{\partial\ln Z_{0}}{\partial R}+\frac{M_{0}^{2}}{2\pi}\big{\{}\gamma_{E}+1+\ln(R/2a)\big{\}}$$

$${\mathcal{E}}_c=-\frac{m_s}{2\pi }{\int }_{-B}^{B}d\theta \cosh (\theta )ϵ(\theta ),$$

$${\int }_{-B}^{B}d\theta {\tilde{K}}_c(\theta -{\theta }^{\prime })ϵ({\theta }^{\prime })=h_1-m_s\cosh (\theta ),$$

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

$$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\Big{(}\frac{\pi\omega}{2}\Big{)}.$$

$${\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)}.$$

$$p_i=\frac{{\alpha }_i^2}{2\pi },i=1,2.$$

Full text

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 QFT-s 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 Feynman (1982); Lloyd (1996); Abrams and Lloyd (1997); Douçot et al. (2004); Büchler et al. (2005); Cirac et al. (2010); Casanova et al. (2011); Jordan et al. (2014, 2012).

Here, we describe analog quantum simulators for some QFT-s with superconducting quantum electronic circuit (QEC) lattices Houck et al. (2012); Devoret and Schoelkopf (2013). Specifically, we are interested in obtaining bosonic QFT-s directly, and not as consequences of mathematical manipulations (bosonization) of an underlying fermionic or spin system Fendley et al. (1995a, b); Milliken et al. (1996); Lesage et al. (1997); Giuliano and Sodano (2005); Gritsev et al. (2007). Even in $1+1$ space-time dimensions, this is crucial to distill the true behavior of a bosonic system from a bosonized fermionic system Bukhvostov and Lipatov (1981); Saleur (1999), despite the well-known fermion-boson correspondences Giamarchi (2003); Gogolin et al. (2004). 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 Fisher et al. (1989); Sachdev (2011). Another scheme is to use the Josephson effect to give rise naturally to the cosine potential Tinkham (2004) 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 Devoret (1997) (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 QFT-s.

The playground for bosonic nonlinear QFT-s 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 with 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 semi-classical description 111This is true also if one tries to “backtrack” to the fermionic formulation, where “non-Fermi liquid” fixed points are then encountered. In fact, one of the simplest new aspects to formulate and study in multi-component bosonic theories is integrability. Integrability leads to a factorized multi-particle scattering matrix Iagolnitzer (1978); Zamolodchikov and Zamolodchikov (1979). This, in turn, allows analytic computation of thermodynamic and transport properties of these QFT-s Zamolodchikov (1990); Smirnov (1992). The ability to do such computations is crucial to understand non-perturbative aspects of strongly interacting systems - for instance, charge fractionalization Senthil and Fisher (2000); Xu and Fu (2010); Roy et al. (2017) and its effect on full-counting statistics Nazarov and Blanter (2009), both topics of high current interest.

It is well known that the bulk and boundary SG models are classically integrable Rajaraman (1982); Ablowitz and Segur (2006). This was largely used historically to establish their quantum integrablity as well  Mandelstam (1975); Zamolodchikov and Zamolodchikov (1979); Ghoshal and Zamolodchikov (1994). 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 Zarembo (2017)- 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 Bernard and LeClair (1991). 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 Belavin et al. (1984); Zamolodchikov (1987, 1989); Francesco et al. (1997); Mussardo (2010).

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

[TABLE]

where $\phi_{1},\phi_{2}$ are two bosonic fields, $M_{0}$ the interaction strength and $\alpha_{1},\alpha_{2}$ 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 Affleck et al. (2001). 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 Ameduri and Efthimiou (1998). 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 Bukhvostov and Lipatov (1981); Fateev (1996); Lesage et al. (1997, 1998); Ameduri et al. (1999); Bazhanov et al. (2018). 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 set-up where such behavior - which, we believe, is the norm rather than the exception in many multi-field 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 Fateev (1996), 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 non-compact 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 The 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 $\varphi_{1},\varphi_{2}$, where $\varphi_{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 Lesage et al. (1998), the two quantum numbers of each soliton scatter independently of each other. The factorized scattering matrix is expected to be given by Fateev (1996) $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$ Mussardo (2010); Zamolodchikov and Zamolodchikov (1979), where $p_{1,2}=\alpha_{1,2}^{2}/2\pi$ Fateev (1996).

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}\varphi_{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 Wiener-Hopf technique Fateev (1996); Fateev et al. (1993) 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 Fateev (1996)), in addition to its usual quadratic dependence:

[TABLE]

Here $m_{s}$ is the mass of the $\varphi_{1}$ soliton sol . It is obtained in terms of the action parameters using conformal perturbation theory (CPT) Klassen and Melzer (1991). We find that $m_{s}=M_{0}/\sin(\pi p_{1}/2)$, which matches the result of Ref. Fateev, 1996. The calculation is is described in detail in Appendix A.

III The double sine-Gordon model at nonzero temperature

The thermodynamics of QFT-s in an infinite volume can be inferred from its scattering matrix Dashen et al. (1969) using the TBA technique Zamolodchikov (1990). 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 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 Zamolodchikov and Zamolodchikov (1979). The TBA analysis is, in general, highly non-trivial, the scattering in each sector being non-diagonal. We consider the simpler case $p_{1}=1/(n+1)$ when the first sector scattering is diagonal. The diagonalization problem of the second sector is done by the Takahashi-Suzuki classification of the solution of the Bethe equations Takahashi and Suzuki (1972); Fendley and Intriligator (1992); Takahashi (2005) and the algebraic Bethe ansatz Korepin et al. (1993); Faddeev (1996). The detailed derivation of the TBA equations is given in Appendix B. The resulting TBA equations have an universal form for all $n$:

[TABLE]

The TBA kernel is $\hat{\varphi}_{n}(\theta)=(n+1)/\cosh[(n+1)\theta]$, while ${\cal 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 Zamolodchikov and Zamolodchikov (1979) and $a\star b=\int a(\theta-\theta^{\prime})b(\theta^{\prime})d\theta^{\prime}/2\pi$. The TBA diagram is given in Fig. 3.

We check our TBA equations by analytic computation of the central charge in the conformal (UV) limit. This yields the expected result of $c_{\rm{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_{\rm{eff}}=-6R^{2}f/\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_{\rm{eff}}$ Zamolodchikov (1991a); Fendley et al. (1993) upon deviations from the conformal limit, in addition to the polynomial corrections in powers of $(m_{s}R)^{2}$ Klassen and Melzer (1991). We calculate these corrections by numerically solving the TBA equations. The logarithmic correction for the different values of $n$ is then verified using CPT. We find the logarithmic correction to $c_{\rm{eff}}$ to be

[TABLE]

More details on the TBA and the explicit forms of the different polynomial corrections to $c_{\rm{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  Zamolodchikov (1991a); Fendley et al. (1993) and should be contrasted with a very different type of logarithmic correction that arises in the Fateev or sinh-Gordon models Fateev (1996); Zamolodchikov (2006).

IV The 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 $\varphi_{1,2}$, the dSG action can be viewed as two SG actions coupled by $B\int d^{2}x\partial_{\mu}\varphi_{1}\partial_{\mu}\varphi_{2}$, where $B=\pi/\alpha_{1}^{2}-\pi/\alpha_{2}^{2}$. Recall that in QEC-s, the bosonic field is time-integral of the voltage at a point. Thus, the first coupling term, $\partial_{t}\varphi_{1}\partial_{t}\varphi_{2}$, is a voltage-voltage coupling, realized in QEC-s with a capacitance ($C$). The second coupling term, $\partial_{x}\varphi_{1}\partial_{x}\varphi_{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 QEC-s for the two SG models are two one-dimensional arrays [see panel (e) of Fig. 2] 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$ Manucharyan et al. (2009); Koch et al. (2009); Le Hur (2012); Goldstein et al. (2013), where $R_{Q}(\simeq 6.5\ 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 are exponentially small $\sim e^{-\sqrt{E_{J}/E_{C_{J}}}}$ Koch et al. (2007). 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/(2/K+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 wavepacket Ablowitz and Segur (2006) 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 wavepacket 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 Guarcello et al. (2018a), 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 panel (e) of Fig. 2 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. Lesage et al., 1997, 1998 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 QFT-s using QEC-s open new possibilities in investigation of their non-perturbative properties. Radio-frequency measurements possible with QEC-s allow a more detailed verification of the QFT-s’ properties. At the same time, QEC simulation of integrable QFT-s may be viewed as a method to benchmark these analog simulators by comparing experiments to analytical theoretical predictions. This leads to a systematic investigation of more general QFT-s, by systematically including perturbations which break integrability. Using the robust controllability and customizability of QEC-s, 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 panel (e) of Fig. 2]. 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 QFT-s may be obtained by specific arrangements of Josephson junctions in tailored geometries Bergeal et al. (2010); Roy et al. (2015); fou . Fermionic modes may potentially also be included using the topological Josephson effect Fu (2010). Finally, the QEC-s proposed in this work can be generalized to simulate higher dimensional QFT-s (integrable and non-integrable) Zamolodchikov (1981).

Acknowledgments

A.R. acknowledges fruitful discussions with Jonathan Belletete, Michel Devoret, Hermann Grabert, Michal Pawelkiewicz, and A. Douglas 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 The 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}\rightarrow\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

[TABLE]

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

[TABLE]

In limit $h_{1}/M_{0},h_{2}/M_{0}\rightarrow\infty$, the interaction term ${\cal S}_{\rm{int}}$ can be neglected. The ground state energy can be computed by completing squares in a Gaussian integral and is given by

[TABLE]

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

[TABLE]

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

[TABLE]

and the ground state energy is

[TABLE]

Next, we look at finite $h_{1}/M_{0},h_{2}/M_{0}$ and consider fluctuations of the fields $\varphi_{1,2}$ around their asymptotic values given by Eq. [8]. This is where our results differ from those obtained in Ref. Fateev, 1996. As will be shown below, the corrections to the limit $h_{1}/M_{0},h_{2}/M_{0}\rightarrow\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

[TABLE]

where $\langle\tilde{\varphi}_{i}\rangle=0$. The original fields $\phi_{1,2}$ can be written has

[TABLE]

where we have defined $\tilde{\phi_{1,2}}=\sqrt{\pi}(\tilde{\varphi}_{1}\pm\tilde{\varphi}_{2})/\alpha_{1}$. Rewriting the action in terms of $\tilde{\phi}_{1,2}$ and considering the case when $Q_{2}=0$, we have ${\cal S}_{0}+{\cal S}_{\rm{mag}}=\tilde{\cal S}_{0}+LR{\cal E}_{c}$, where

[TABLE]

Here, $L$ is the spatial extent of the 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 ${\cal S}_{\rm{int}}$ perturbatively. We consider the case when $Q_{2}=0$. The interaction term can be written as

[TABLE]

This leads to

[TABLE]

where we have used the fact that $\langle{\cal S}_{\rm{int}}\rangle_{0}=0$, $Z_{0}=\int{\cal D}\hat{\phi}_{1}{\cal D}\hat{\phi}_{2}{\rm{exp}}(-\hat{\cal 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 Zamolodchikov (1995)

[TABLE]

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

[TABLE]

where $\gamma_{E}$ is the Euler constant and $a$ is a lattice cut-off. Using this result, the specific ground state energy, up to second order in $M_{0}$, is given by

[TABLE]

where

[TABLE]

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 $\varphi_{1}$ fields acquire an additional energy $-h_{1}Q_{1}$ and as long as $h_{1}$ is larger than the mass of the $\varphi_{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 wavevector’. The specific ground state energy is then given by Polyakov and Wiegmann (1983); Wiegmann (1985); Korepin et al. (1993)

[TABLE]

where $m_{s}$ is the mass of the $\varphi_{1}$ soliton and the quasiparticle energies $\epsilon(\theta)$ satisfy

[TABLE]

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

[TABLE]

where $\delta_{p_{i}}$-s are soliton-soliton scattering phase-shifts (for instance, see Ref. Fateev, 1996 for explicit forms). The Fourier transform of $\tilde{K}_{c}$ is given by

[TABLE]

The leading order energy contribution in the limit of $h_{1}/M_{0},h_{2}/M_{0}\rightarrow\infty$ can be obtained by the standard manipulations Fateev (1996); Fateev et al. (1993) and the result is identical to that of the Fateev model. This leads to

[TABLE]

Comparing Eqs. [10, 22], we get

[TABLE]

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}\rightarrow\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

[TABLE]

The kernel factorizes as $K_{c}(\omega)=1/N(\omega)N(-\omega)$ Fateev (1996), where

[TABLE]

where

[TABLE]

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

[TABLE]

where $m_{s}$ is the mass of physical excitations of solitons of $\varphi_{1}$, $\alpha(\omega)=N(\omega)/N(-\omega)$ and

[TABLE]

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

[TABLE]

where $C^{\prime}_{+}$ 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

[TABLE]

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

[TABLE]

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

[TABLE]

The contributions to the ground state energy can be identified as follows. The first term on the right hand side is ${\cal E}_{c}$, the leading order contribution to the ground state energy when $h_{1}/M_{0},h_{2}/M_{0}\rightarrow\infty$. In the second term, the $h_{1}$-independent term can be identified as the ground state energy ${\cal 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

[TABLE]

where

[TABLE]

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

[TABLE]

Furthermore, comparing the coefficient of $\ln h_{1}$ term in Eqs. (16, 31) leads to an exact expression connecting the interaction parameter $M_{0}$ of the action with the masses of the $\varphi_{1}$ solitons:

[TABLE]

Here, we have used the definition of $Q_{1}$ in Eq. [9] and the relation Eq. [23]. 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 The 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}\geq 1$ Fateev (1996). 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,\ldots,n$ Zamolodchikov and Zamolodchikov (1979). The TBA analysis of this problem is in general highly non-trivial, being each scattering sector source of non diagonally. 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 consists 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 Takahashi and Suzuki (1972); Fendley and Intriligator (1992); Takahashi (2005).

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:

[TABLE]

where $\ i=1,\ldots,n$. In this equation, $\theta$ is the rapidity of the particle going around the circle, $\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 Zamolodchikov and Zamolodchikov (1979), 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 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 will be needed in the identification of the role of the massive physical particles. This is done using the algebraic Bethe ansatz technique Korepin et al. (1993); Fendley and Intriligator (1992); Faddeev (1996). 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 manuscript (Appendix D) and here, we present directly the final results:

[TABLE]

In the above equation, the TBA kernel is given by $\hat{\varphi}_{n}(\theta)=(n+1)/\cosh[(n+1)\theta]$, ${\cal 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^{\prime}}{2\pi}a(\theta-\theta^{\prime})b(\theta^{\prime})$. The pictorial representation for the TBA equations is given Fig. 3. 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. [B], 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 ${\cal L}$-s and $L$-s, which allows us to write the TBA equations concisely in terms of a single universal kernel $\hat{\varphi}_{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. Fateev, 1996, 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

[TABLE]

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\rightarrow\infty$) and IR ($T\rightarrow 0$) limits Zamolodchikov (1991a, b). The resultant expression for the central charge can be written as

[TABLE]

where ${\cal L}_{\rm{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,\ldots,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 and so, they do not contribute to the expression of the central charge Zamolodchikov (1991a). The solutions of the TBA equations in the UV limit are given by

[TABLE]

while for the IR limit, the non-zero $y$-s are given by

[TABLE]

Plugging these solutions in Eq. [37] and using some remarkable dilogarithm identities Lewin (1981), we arrive the desired result of

[TABLE]

Next, we compute the effective central charge as we move away from the UV limit. The perturbation ${\cal S}_{\rm{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 Zamolodchikov (1991a); Fendley et al. (1993) as one moves away from the UV limit, in addition to the polynomial corrections in powers of $(m_{s}R)^{2}$ Klassen and Melzer (1991). 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

[TABLE]

where the kernels are given by

[TABLE]

We solve Eq. [B] numerically and plug it in Eq. [36] to obtain the effective central charge $c_{\rm{eff}}$. The result is shown in Fig. 4.

The first six correction terms capturing the deviation from the UV limit is obtained by fitting the numerical data to the following expression:

[TABLE]

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

[TABLE]

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

[TABLE]

The coefficient of the logarithmic term in the effective central charge calculation can be checked analytically using CPT. Taking Eq. [A] for $h_{1}=0$ and using Eq. [33], we arrive at the expression for the free energy:

[TABLE]

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

[TABLE]

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. [B, B]). We have checked our results for $n=4$ and $n=5$ as well.

Appendix C The double sine-Gordon model with quantum circuits

In this section of the 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. 5 (a). Define the node-flux at a point in space-time $(x,t)$ as $\phi(x,t)$ Devoret (1997), which plays the role of the bosonic field. As shown in Refs. Clerk et al., 2010; Goldstein et al., 2013; Roy and Devoret, 2016, 2018, 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 ${\cal S}_{\rm{LL}}=\frac{1}{2\pi K}\int d^{2}x(\partial_{\mu}\phi)^{2}$, where we have rescaled the space and time axis to set the plasmon velocity to be 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 ${\cal S}_{\rm{SG}}={\cal S}_{\rm{LL}}-E_{J_{0}}\int d^{2}x\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\leq(\geq)\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. 5 (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 ${\cal S}_{{\rm{SG}},i}$, $i=1,2$:

[TABLE]

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

[TABLE]

while the second leads to

[TABLE]

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 ${\cal S}_{{\rm{SG}},1}+{\cal S}_{{\rm{SG}},2}$ $+{\cal S}_{M}+{\cal S}_{C}$. From this action, simple algebraic manipulations give rise to the dSG action in terms of the fields $\varphi_{1,2}$, given by

[TABLE]

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,

[TABLE]

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

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 finite number of such unit cells. There are two main difficulties in realizing arrays 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}}}}$ Manucharyan et al. (2009); Goldstein et al. (2013). The same problem also provides a limitations on the frequencies over width 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}}$, the plasmon velocity Goldstein et al. (2013). 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}$ Manucharyan et al. (2009). Despite these difficulties, arrays with up to ${\cal O}(10^{5})$ Josephson junctions have been reliably built and experimentally analyzed Kuzmin et al. (2018) 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 Fominaya et al. (1997). 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 Fominaya et al. (1997); Rabani et al. (2008). A detailed theoretical analysis for the specific heat measurements for the sine-Gordon model has been made in Refs.  Guarcello et al., 2018a, b.

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 Hida and Eckern (1984); Hida (1985) indicate a finite width for the sine-Gordon minima which presumably leads to a finite width for the different particles in the spectrum. 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

[TABLE]

where $\theta$ is the rapidity of the particle going around the world, $\theta_{k_{s,a,b_{1},b_{2}}}$ the rapidities of the external “incoming” solitons, antisolitons and breathers and in which $\lambda^{p_{2}}_{i}$ 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^{p_{2}}_{SG}$. The diagonalization procedure for repulsive sine-Gordon 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 equations (LABEL:periodicity) the only care will be needed in the identification of the role of the massive physical particles. Through the algebraic Bethe ansatz technique Zamolodchikov (1991a); Fendley and Intriligator (1992), the eigenvalues of the repulsive sine-Gordon transfer matrix are

[TABLE]

The rapidities $y_{r}$ label the Bethe eigenvectors of the transfer matrix and have to satisfy the Bethe equation

[TABLE]

The presence of the external particles does not invalidate the string classification Takahashi and Suzuki (1972); Takahashi (2005) 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 equation (53) is equivalent to the classification of the $y^{\prime}_{r}=(1-p_{2})y_{r}$ of an equation like (53) but with $p^{\prime}_{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)}$

[TABLE]

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

[TABLE]

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 (54) in the Bethe equation (53) and by taking the imaginary part of the logarithmic derivative with respect to the external rapidity. The result for the pseudoparticles is

[TABLE]

where

[TABLE]

are the kernels of the particle $i$ with the physical particle denoted with [math] and with the other pseudoparticles $j$. The convolution normalization is $a\star b=\int\frac{d\theta^{\prime}}{2\pi}a(\theta-\theta^{\prime})b(\theta^{\prime})$. For sine-Gordon 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 the equation (56) to the equation

[TABLE]

where the breather kernels are

[TABLE]

The equations for the densities of the physical particles states can be obtained by taking the logarithmic derivative of equations (LABEL:periodicity) with the explicit formulas for the eigenvalues (52) and (54). The result is

[TABLE]

The new kernels are

[TABLE]

where we denote with

[TABLE]

the soliton-soliton kernel of the pure $p_{2}$-sineGordon 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 dke^{ik\theta}\tilde{\phi}(k)$ and assume the form

[TABLE]

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 results into a drastic simplification of the system

[TABLE]

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

[TABLE]

and the entropy, in the Stirling approximation,

[TABLE]

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

[TABLE]

Taking into account the constraints of Eq. (65), we end up to the following TBA system for the pseudoenergies $\epsilon_{i}(\theta)$ defined as $\rho^{+}_{i}/P_{i}=e^{-\epsilon_{i}}/(1+e^{-\epsilon_{i}})$

[TABLE]

where $R$ is the inverse temperature. Here we have introduced the notation $L_{i}(\theta)=\log(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

[TABLE]

The functions ${\cal L}_{i}$ are the positive part of $\epsilon_{i}$, i.e. ${\cal L}_{i}=\log(1+e^{\epsilon_{i}})$. The universal kernel is $\tilde{\varphi}_{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. (B).