# Quantum simulation of the tricritical Ising model in tunable Josephson
junction ladders

Lorenzo Maffi  Center for Quantum Devices and Niels Bohr International
Academy, Niels Bohr Institute, University of Copenhagen, DK–2100 Copenhagen,
Denmark  Niklas Tausendpfund  Forschungszentrum Jülich GmbH, Institute of
Quantum Control, Peter Grünberg Institut (PGI-8), 52425 Jülich, Germany
Institute for Theoretical Physics, University of Cologne, D-50937 Köln,
Germany  Matteo Rizzi  Forschungszentrum Jülich GmbH, Institute of Quantum
Control, Peter Grünberg Institut (PGI-8), 52425 Jülich, Germany  Institute for
Theoretical Physics, University of Cologne, D-50937 Köln, Germany  Michele
Burrello  Center for Quantum Devices and Niels Bohr International Academy,
Niels Bohr Institute, University of Copenhagen, DK–2100 Copenhagen, Denmark

###### Abstract

Modern hybrid superconductor-semiconductor Josephson junction arrays are a
promising platform for analog quantum simulations. Their controllable and non-
sinusoidal energy/phase relation opens the path to implement nontrivial
interactions and study the emergence of exotic quantum phase transitions.
Here, we propose the analysis of an array of hybrid Josephson junctions
defining a two-leg ladder geometry for the quantum simulation of the
tricritical Ising phase transition. This transition provides the paradigmatic
example of minimal conformal models beyond Ising criticality and its
excitations are intimately related with Fibonacci non-Abelian anyons and
topological order in two dimensions. We study this superconducting system and
its thermodynamic phases based on bosonization and matrix-product-states
techniques. Its effective continuous description in terms of a three-frequency
sine-Gordon quantum field theory suggests the presence of the targeted
tricritical point and the analysis of order parameters and correlation lengths
confirm this picture. Our results indicate which experimental observables can
be adopted in realistic devices to probe the physics and the phase transitions
of the model. Additionally, our proposal provides a useful one-dimensional
building block to design two-dimensional scalable Josephson junction arrays
with exotic topological order.

The rapid advances in the fabrication of superconducting/semiconducting
heterostructures [1, 2] allow for the realization of Josephson junction arrays
(JJAs) with an unprecedented control and tunability over their physical
parameters [3, 4, 5]. State-of-the-art electron beam lithography and etching
techniques enable the realization of superconducting (SC) arrays with
exquisite geometrical precision and scalability. Epitaxial growth consents to
create pristine interfaces between a semiconducting substrate and SC islands,
thus providing the possibility of controlling these setups through voltage
gates. These fabrication developments are flanked by remarkable advances in
the measurement techniques which include microwave spectroscopy to study the
1D strongly correlated systems emerging in Josephson junction chains [6, 7]
and transport measurements to investigate the intricate thermodynamic
properties of these systems [8, 3, 4, 5, 7]. Such progresses have brought JJAs
right back into the arena of analog quantum simulation platforms, where they
started their journey decades ago. The simultaneous tunability of the junction
transparencies [2, 9, 10, 11] and magnetic fluxes opens indeed the path to
tailor models of interest, among which quantum field theories and integrable
models [12, 13, 14]. In particular, the experimental achievement of
multicritical points, with peculiar conformal field theories (CFTs) associated
to them [15], becomes at reach [16].

In this work we formulate a blueprint for the quantum simulation of the
tricritical Ising (TCI) CFT in a tunable Josephson junction ladder. The
reasons of interest for this model are multiple. It constitutes the simplest
example of CFT beyond the Ising model, and its particle content includes
excitations that share the same fusion properties of Fibonacci non-Abelian
anyons. Successfully implementing this model will open the way to engineer
exotic topological order in 2D arrays in the spirit of the wire constructions
of Refs. [17, 18, 19, 20]. Moreover, the TCI model stands as a strong
potential candidate to observe the emergence of supersymmetry [21, 22, 23].
Notably, to our knowledge, no experimental realization of a quantum
tricritical Ising phase transition in 1D has ever been observed, nor have its
critical exponents been measured.

Indeed, the quantum simulations of CFTs beyond the Ising universality class
face both experimental and theoretical challenges: the most recent theoretical
proposals rely on advanced constructions based on Majorana modes [22, 24, 25,
23, 26, 20], extended Hubbard models with staggering potentials [27, 28] or
nontrivial mappings between microscopic lattice operators and the field
content of the CFTs [29]. In this context, the main mechanism to achieve a TCI
point is to consider physical platforms like Rydberg atom systems [30, 31] and
ultracold atoms in tilted optical superlattices [32] that are described by
discrete models with a continuous Ising phase transition turning into a first-
order phase transition (FOPT) at the tricritical point.

A significant advancement offered by JJAs is the ability to provide a feasible
way to directly implement nontrivial interacting bosonic quantum field
theories (QFTs) [13, 16]. We make use of this aspect to present a ladder
system that embodies a three-frequency sine-Gordon model and, as we will show,
can be tuned to naturally flow towards the TCI point at low energy. The chosen
ladder geometry offers an alternative construction compared to previous works
on superconducting chains [14, 16], and opens a path to the realization of 2D
devices with non-Abelian topological order. To achieve our goal, we utilize a
blend of analytical techniques, including mean field analysis and bosonization
[33], complemented by numerical analysis based on variational uniform matrix
product state ansatz (VUMPS) [34, 35, 36].

#### The triple Josephson junction.-

The building block of our 1D construction consists of two E-shaped SC islands
facing each other and grown on top of a semiconducting substrate [Fig. 1(a)].
Schematically, we can model this element as three parallel Josephson junctions
(JJs) where Andreev bound states induced in the semiconductor mediate the
Cooper pair tunneling [37, 38]. For simplicity we assume that each junction is
defined by a single transport channel with transparency $T_{p}$ ($p=1,2,3$)
(the case of multichannel junctions is briefly discussed in [39]). In the
short junction limit, the energy/phase relation of each single-channel JJ $p$
can be written as [37]

$\mathcal{E}_{J}^{(p)}\left(\varphi\right)=-\Delta\sqrt{1-T_{p}\sin^{2}{\left(\varphi/2\right)}},$
(1)

where $\varphi$ is the phase difference between the two islands and $\Delta$
is the SC gap induced by proximity in the semiconducting substrate. We assume
$\Delta$ to be the largest energy scale in our system. According to equation
(1), high-transparency channels in the semiconductor favour the onset of the
coherent tunneling of multiple Cooper pairs [40]. These events correspond to
the higher harmonics contribution, $\cos(n\varphi)$ with $n>1$, to the
Josephson energy/phase dispersion. The transparencies of the external
junctions, $T_{1}$ and $T_{3}$, are controlled by external electrostatic gates
[Fig. 1(a)] which can vary the carrier density in the surrounding
semiconductor.

Figure 1: (a) Two E-shaped SC islands are connected through three parallel
junctions characterized by two magnetic fluxes $\Phi_{1}/2\pi$ and
$\Phi_{2}/2\pi$ in units of $\Phi_{0}=hc/2e$. The external junctions are
controlled by two electrostatic gates at potential $V_{G1}$, $V_{G3}$. This
triple JJ element allows us to generate and control the local potential (3) at
each rung of the ladder geometry (b). We connect the SC islands along the legs
through a standard sinusoidal JJ. The fluxes of the triple JJ elements are
staggered along the ladder [39]. Mutual rung capacitances and the island self-
capacitances determine the electrostatic interactions $V_{\perp}$ and $E_{C}$.

Additionally, the geometry of the triple JJs defines two loops in which
controllable magnetic fluxes $\Phi_{1}/2\pi,$ and $\Phi_{2}/2\pi$ can be
inserted [Fig. 1(a)]. Therefore, Cooper pairs tunneling across different paths
acquire relative Aharanov-Bohm phases leading to interference. With a suitable
gauge choice, we write the full Josephson energy of one triple JJs element as
the sum of three contributions [39]:

$V_{J}\left(\varphi\right)=\leavevmode\nobreak\
\mathcal{E}_{J}^{(1)}\left(\varphi-\Phi_{1}\right)+\leavevmode\nobreak\
\mathcal{E}_{J}^{(2)}\left(\varphi\right)+\mathcal{E}_{J}^{(3)}\left(\varphi+\Phi_{2}\right)\
.$ (2)

In the following, we set the same flux inside the two plaquettes
$\Phi_{1}=\Phi_{2}=\Phi$, and the same transparencies for the two external
junctions $T_{1}=T_{3}$ (i.e., the same gate voltages $V_{G1}$ = $V_{G3}$ for
the case of single-channel junctions). Thus, the exchange of the SC islands
(or, equivalently, $\varphi\to-\varphi$) corresponds to a
$\mathds{Z}_{2}$-symmetry, which is necessary for the multicritical Ising
physics to emerge. Such symmetry is embodied in a measurable odd current/phase
relation of the triple JJ. The presence of multiple channels in the junctions
or unequal plaquette areas may, however, cause perturbations that explicitly
break this $\mathbb{Z}_{2}$ symmetry [39].

These perturbations would prevent observing the critical features of the model
whenever the energy gaps they open are larger than the energy resolution
achievable in experiments, which is given by the finite size $L$ and the
finite temperature $T$.

In the symmetric setup, the total Josephson potential can be expanded as

$V_{J}\left(\varphi\right)=\mu_{1}({\bf X})\cos{\left(\varphi\right)}+\\\
\mu_{2}({\bf X})\cos{\left(2\varphi\right)}+\mu_{3}({\bf
X})\cos{\left(3\varphi\right)}+\ldots.$ (3)

The Fourier coefficients $\mu_{n}$ [39] depend on the values of the external
parameters ${\bf
X}=\left(T_{1}\cos{\left(\Phi\right)},\,T_{1}\sin{\left(\Phi\right)},\,T_{2}\right)$
which span a solid cylinder.

In the following, we will adopt many copies of this triple JJ to build a 1D
ladder geometry, thus promoting the phase difference $\varphi$ to a position-
dependent field. In light of this, a preliminary mean-field analysis allows us
to gain a qualitative understanding on the onset of a TCI point by
investigating the Josephson potential $V_{J}\left(\varphi\right)$ as function
of ${\bf X}$. In a semiclassical picture, a tricritical point can indeed be
obtained when three potential minima merge into one [41, 42, 43]. In the
potential landscape defined by $V_{J}(\varphi)$ with
$\varphi\in\left(-\pi,\pi\right]$, for any section of the cylinder at fixed
$T_{2}$, there exists a point $\left(T_{1},\Phi\right)_{c}$ where this merging
happens and $V_{J}(\varphi)$ is approximated by a $\varphi^{6}$ local
potential. This fact suggests the first connection to the tricritical Ising
model and its Ginzburg-Landau (GL) formulation [41, 42, 43].

Figure 2: Given $\varphi_{\rm min}$ the global minimum of the Josephson
potential (3), we report the value of $\left|\sin{\left(\varphi_{\rm
min}\right)}\right|$ for each point of the parameter space with fixed
$T_{2}=0.6$. This function plays the role of a mean-field magnetization and is
odd under the $\mathds{Z}_{2}$-symmetry. The first harmonic in Eq. (3) mostly
prevails and picks $\varphi_{\rm min}=\left\\{0,\pi\right\\}$ depending on its
sign (regions I and III respectively). These regions correspond to
$\mathds{Z}_{2}$-symmetric configurations. The competition with higher
harmonics gives rise to an ordered phase with two degenerate global minima
$\varphi_{\rm min}=\pm\varphi_{0}$ (region II). In the zoomed area showed in
the inset, the transition between region I and II can be either discontinuous
with three degenerate minima (yellow line) or continuous with the merging of
the two minima in $\varphi_{\rm min}=0$. At the meeting point we have a
tricritical point around which a three well potential
$V_{J}=g_{2}\varphi^{2}+g_{4}\varphi^{4}+\varphi^{6}$ approximates Eq.(3). The
dashed black line indicates the points where $g_{4}=0$ and new local minima
emerges.

#### 1D model.-

We design a 1D quantum simulator to achieve a TCI point by arranging a set of
identical triple JJs with potential $V_{J}$ in parallel, as depicted in Fig.
1(b), in order to implement a multi-frequency sine-Gordon model at low
energies. In particular, we consider a regime with three sine-Gordon
interactions [44], stemming from the potential (3). The Hamiltonian of the JJ
ladder [Fig. 1(b)] is:

$\begin{split}\widehat{H}=\sum_{j=0}^{L-1}\Biggl{[}\sum_{\alpha={a},{b}}&\\!\left(E_{C}\widehat{N}_{\alpha,j}^{2}-E_{J}\cos{\left(\hat{\varphi}_{\alpha,j+1}-\hat{\varphi}_{\alpha,j}\right)}\right)\Biggr{.}\\\
&+\Biggl{.}V_{\perp}\,\widehat{N}_{{a},j}\widehat{N}_{{b},j}+V_{J}\left(\hat{\varphi}_{{a},j}-\hat{\varphi}_{{b},j}\right)\Biggr{]},\end{split}$
(4)

where $\hat{\varphi}_{\alpha,j}$ represents the phase operator of the $j$-th
island on the leg $\alpha\in\left\\{{a},{b}\right\\}$. Along the legs, the SC
islands are connected through JJs in a standard sinusuoidal regime with
Josephson energy $E_{J}$. This energy scale can vary from $E_{J}\simeq
h50\;$GHz [10] down to $E_{J}=0$ for completely depleted junctions. The
dynamics of the SC phases in Eq. (4) is dictated by charging effects,
described by the charge operators $\widehat{N}_{\alpha,j}$, canonically
conjugated to the SC phases,
$[\widehat{N}_{\alpha,j},e^{i\hat{\varphi}_{\alpha,j}}]=-e^{i\hat{\varphi}_{\alpha,j}}$.
We consider in particular an on-site electrostatic repulsion $E_{C}$ and a
rung repulsive interaction $V_{\perp}$.

To obtain the rung potentials $V_{J}$ in Eq. (4), the pattern of magnetic
fluxes in the system must be carefully considered: a uniform magnetic field
breaks time-reversal invariance driving the system into Meissner chiral phases
[45, 46, 47, 48, 49, 50, 51] and does not fulfill the $\mathbb{Z}_{2}$
symmetry on each rung. We consider instead staggered fluxes $(-1)^{j}\Phi$
alternating at each triple JJ [Fig. 1 (b)]. This choice yields the local
effective potential (3) and avoids additional fluxes between subsequent rungs
[39].

The aimed description in terms of the multi-frequency sine-Gordon model
emerges in the regime in which the rung potentials $V_{J}$ and the Josephson
energy $E_{J}$ dominate over the charging effects $E_{C}$ and $V_{\perp}$. In
this Josephson-dominated regime, the system lies away from Mott insulating
phases [52, 53, 48] and phase delocalization due to charge disorder [54, 55,
56] is strongly irrelevant. In this limit, the low-energy physics of the
Cooper pairs in long ladders can be described through bosonization [33]. We
introduce the quantum fields $\hat{\theta}_{\alpha}(x)$ that represent the
island charge density fluctuations
$\widehat{N}_{\alpha,j}/a\approx-\partial_{x}\hat{\theta}_{\alpha}(x)/\pi$,
with $x=ja$ and $a$ the ladder lattice spacing. These fields are dual to the
compact phase fields defined as
$\hat{\varphi}_{\alpha}(x)\approx\hat{\varphi}_{\alpha,j}$.

By introducing the customary charge $c$ and spin $s$ sectors,
$\hat{\varphi}_{c/s}(x)=\left(\hat{\varphi}_{{a}}(x)\pm\hat{\varphi}_{{b}}(x)\right)/\sqrt{2}$,
the Hamiltonian (4) is approximated by:

$\widehat{H}=\sum_{q=c,s}u_{q}\int\dfrac{dx}{2\pi}\left[K_{q}\left(\partial_{x}\hat{\varphi}_{q}\right)^{2}+\dfrac{1}{K_{q}}\left(\partial_{x}\hat{\theta}_{q}\right)^{2}\right]\\\
+\int\frac{dx}{a}\,\sum_{n=1}^{3}\mu_{n}\cos{\left(\sqrt{2}n\hat{\varphi}_{s}\right)}.$
(5)

Eq. (5) describes the two sectors of the model as Luttinger liquids (LLs),
with Luttinger parameters $K_{c/s}\approx\pi\sqrt{E_{J}/\left(2E_{C}\pm
V_{\perp}\right)}$ [48, 39]. The rung potential $V_{J}$ affects only the spin
sector and determines the onset of the targeted multiple sine-Gordon
interactions. To this purpose, the three potential terms in Eq. (5) must be
relevant in the renormalization group sense and induce order in the SC phase
difference $\hat{\varphi}_{s}$, driving the spin sector away from the LL
phase. This sets the constraint $K_{s}>9/4$, which, indeed, is fulfilled for
sufficiently large Josephson energies, when the semiclassical description of
the field $\hat{\varphi}_{s}$ is most accurate. Higher frequency harmonics,
instead, are less relevant and characterized by smaller amplitudes, such that
we neglect them [39].

The interplay of the three sine-Gordon terms $\mu_{n}$ opens the path towards
nontrivial phase transitions [57, 44, 16] between the low-energy massive
phases of the spin sector. In particular, an Ising critical line meets a FOPT
in a tricritical point characterized by the TCI CFT with central charge
$c=7/10$ [44, 16]. Despite the absence of a non-perturbative mapping between
the bosonic fields and the operator content of the TCI, the
$\mathds{Z}_{2}$-symmetry $\hat{\varphi}_{s}\to-\hat{\varphi}_{s}$ allows us
to define convenient observables. In the following, we benchmark the field
theory predictions with numerical tensor-network simulations of the ladder
Hamiltonian (4).

Figure 3: (a): The non-zero expectation value of the single-particle current
$\hat{J}^{(2e)}_{\perp}$ indicates the onset of the ordered phase II. Green
stars mark a discontinuous change of the log-fidelity per site (Eq. (6))
denoting the FOPT between the disordered phases I and III. Panel (b) displays
the discontinuity of $\mathcal{F}$ between the ordered phase II and the
disordered phase I at $X_{2}=0.52$. The order operator
$\hat{J}^{(2e)}_{\perp}$ shows as well a discontinuous jump in correspondence
of the cut b). (c): along the cut c) at $X_{2}=0$, the fidelity susceptibility
$\chi_{\mathcal{F}}$ exhibits a singular behaviour, and the order parameter
changes smoothly, both indicating the second order nature of the transition.
Panel (d) shows the collapse of the correlation length $\xi_{s}$ at $X_{2}=0$
for five values of the bond dimension $D$ by employing the finite entanglement
scaling [58]: $\xi_{s}=D^{\kappa}f\left((X_{1}-X_{1c})D^{\kappa/\nu}\right)$
[39]. This method identifies the critical point at $X_{1c}\simeq-0.288$ and
provides the critical exponents $\kappa=0.79(6)$ and $\nu=1.0(1)$. (e): for
several horizontal cuts with $0.42<X_{2}<0.49$, the critical exponent $\beta$
obtained by fitting the scaling of $\hat{J}^{(2e)}_{\perp}$ exhibits the
emergence of two plateaux close to the Ising and TCI predictions, $\beta_{\rm
IS}=1/8$ and $\beta_{\rm TCI}=1/24$. The bond dimension is fixed at $D=600$
and the fitting procedure does not take into account any finite entanglement
effect.

#### Observables and results.-

We study the phase diagram of our model by using the variational uniform
matrix product state ansatz (VUMPS), [34, 35, 36], to find the ground-state of
the Hamiltonian (4) in the thermodynamic limit. The VUMPS is based on a two-
site elementary cell representing two SC islands on the same rung. The local
Hilbert space is constructed from the charge basis defined by
$\widehat{N}_{\alpha={a}/{b},j}$. For numerical purpose, we truncate its basis
by introducing a cutoff, $\left|N_{\alpha,j}\right|<N_{\rm max}$, with $N_{\rm
max}\geq 6$. We checked the quality of this approximation by monitoring the
occupancy probability of higher charge states [39]. In this basis, the vertex
operators $e^{\pm i\hat{\varphi}_{\alpha,j}}$ in Eq.(4) are represented as
raising/lowering operators.

We set $E_{C}/E_{J}=0.4$ and $V_{\perp}/E_{J}=0.65$, corresponding to
$K_{s}\approx 8$. This favours the clean emergence of the transition lines as
the interactions are strongly relevant, yielding sizeable energy gaps in the
spin sector. The Fourier components $\mu_{n}$ in Eq. (3) are determined from
Eq. (1) with an induced gap of $\Delta/E_{J}=50$ and $T_{2}=0.6$, consistent
with Fig. 2 [39].

We identify the phases of the model with the labels I, II and III as in Fig.
2, and, in order to distinguish them, we employ the local order operator
$\hat{J}^{(2e)}_{\perp}(x)=\sin{\left(\sqrt{2}\hat{\varphi}_{s}(x)\right)}$
representing the normalized single-particle contribution to the total rung
current in units of $\mu_{1}$. In the VUMPS simulations, the symmetry-broken
phase II is signalled by a finite $\braket{\hat{J}^{(2e)}_{\perp}}$ [Fig.
3(a)], and it aligns with the mean-field predictions in Fig. 2. The symmetric
phases I and III broaden away from the semiclassical limit due to the dominant
scaling behavior of the first-harmonic interaction. The order parameter allows
us, in particular, to investigate the boundary between the disordered phase I
and the ordered phase II: a neat jump in $\braket{\hat{J}^{(2e)}_{\perp}}$
marks a FOPT for $X_{2}=T_{1}\sin{\left(\Phi\right)}\gtrsim 0.475$, as
exemplified in Fig. 3(c), while a gradual change in the region
$\left|X_{2}\right|\lesssim 0.475$ indicates the onset of a continuous
transition, as exemplified for $X_{2}=0$ in Fig. 3(d). For
$\left|X_{2}\right|\gtrsim 0.55$, instead, the ordered phase closes, and we
observe a FOPT between phases I and III, consistently with the semiclassical
predictions.

This picture is confirmed by the analysis of the ground state fidelities [59,
60, 61, 62]. Given the abrupt change of the ground state $\ket{\psi\left({\bf
X}\right)}$ across the FOPT, the average log-fidelity per site [61]

$\mathcal{F}\left(\bf
X,\delta\right)=-\lim_{N\to\infty}\dfrac{1}{N}\log\left(\braket{\psi({\bf
X}-\delta)}{\psi({\bf X}+\delta)}\right),$ (6)

displays a clean discontinuity [Fig. 3(c)], at fixed $\delta$. On the other
hand, across the lower cut the fidelity susceptibility
$\chi_{\mathcal{F}}=\mathcal{F}/\delta^{2}$ shows a more gradual singular
behaviour and exhibits the typical peak of a second-order phase transition in
Fig. 3(d) [39].

The universal collapse of the spin correlation length $\xi_{s}$ according to
finite entanglement scaling ansatz [58, 39] gives strong indication of an
Ising phase transition at $X_{2}=0$ [Fig. 3(e)]. This method allows us to
precisely locate the critical point $X_{1c}$ and to extrapolate the infinite
bond dimension estimate of the critical exponent $\nu=1.0(1)$ in agreement
with the Ising CFT prediction $\nu_{\rm IS}=1$.

Additionally, our analysis reveals the scaling of the effective magnetization
$\braket{\hat{J}^{(2e)}_{\perp}}\leavevmode\nobreak\ \sim\leavevmode\nobreak\
\left|X_{1}-X_{1c}\right|^{\beta}$, with the critical exponent
$\beta=0.125(3)$ [39] along the $X_{2}=0$ line, compatible with the Ising
prediction $\beta_{\rm IS}=1/8$.

These findings strongly support the classification of the phase transition at
$X_{2}=0$ within the Ising universality class.

By increasing $X_{2}$, the Ising critical line leaves place to the first-order
discontinuity, suggesting the correctness of the semiclassical prediction and
the presence of a TCI point joining the two kinds of phase transitions.

To investigate the onset of the TCI point, we extract the critical exponent
$\beta$ for several horizontal cuts with $X_{2}<0.52$ [below the cut b) in
Fig. 3(a)]. Notably, we observe a good agreement with the Ising prediction for
lower cuts up to $X_{2}\simeq 0.43$. For larger $X_{2}$, the critical exponent
decreases and exhibits a plateau close to the expected TCI value $\beta_{\rm
TCI}=1/24$. Further increasing $X_{2}$ results in a vanishing $\beta$, as
expected by the neat jump of the order parameter corresponding to a FOPT. We
consider this result as a hint of the emergence of the TCI CFT in the phase
diagram of the ladder.

In our simulations finite entanglement effects are accentuated by the massless
LL in the charge sector with $c=1$ throughout the entire phase diagram. These
factors have a quantitative impact on our results, especially making the
central charge estimates along the second-order transition line very noisy
[39]. The measured central charges in correspondence of the plateau observed
in Fig.3(e) are compatible with the presence of a tricritical point.

#### Experimental observables.-

In experimental setups, transport features can be used to explore the phase
diagram of the model. First of all, the central charge of the model can be
directly measured based on heat transport [63]. The heat current across out-
of-equilibrium 1D systems at criticality is known to be proportional to the
central charge of the related CFT at low temperatures [64, 65]: $q=\frac{\pi
ck_{B}^{2}}{12}\left(T_{L}^{2}-T_{R}^{2}\right)$, where $T_{L/R}$ are the
temperatures of two external reservoirs. In our model, this implies that the
heat transport of the symmetric and symmetry-broken phases are entirely
determined by the charge sector of the system ($c=1$). The second-order Ising
line, instead, is characterized by an additional contribution from the spin
sector at criticality, resulting in $c=1+1/2$. This contribution increases to
$c=1+7/10$ at the tricritical point and vanishes on the FOPT. In realistic
devices, finite size and temperature effects determine the profile of the heat
conductance as a function of the system parameters. Nevertheless, a non-
monotonic behavior of the heat conductance across the second-order phase
transition lines and in proximity of the TCI point would provide strong
evidence of the emergence of the related CFTs (see Ref. [66, 27] for the
related study of the central charge).

Another important feature that can be experimentally probed is offered by the
noise in the rung currents. These currents have vanishing expectation values
in all phases, but display long-range order in the symmetry-broken phase.
Accordingly, the power spectrum of their noise constitutes a tool to
distinguish ordered and disordered phases, and its behavior at criticality
provides information on the critical exponents of the model. Finally,
microwave spectroscopy of JJAs [6, 7] allows for the study of the excitation
spectra of the system and can be used to verify the predictions of the TCI CFT
spectra [67, 43, 68, 69, 70]

#### Conclusions.-

We designed a JJ ladder to realize a quantum simulator for the tricritical
Ising CFT. Our construction is based on the properties of hybrid
semiconducting-superconducting Josephson junctions and their non-sinusuoidal
energy/phase relation. In particular, we engineered a triple JJ that allows us
to tune the higher harmonics characterizing this dispersion and we adopted
them to realize a quantum simulator for a multi-frequency sine-Gordon quantum
field theory [44].

The connection between our JJ ladder and this theory has been verified through
bosonization and tensor-networks simulations. Our analysis of the related
phase diagram showed the presence of an ordered phase and demonstrated the
existence of a critical Ising plane connected to a first-order transition in
correspondence of a tricritical Ising line in a three-parameter space.

Our construction does not require the introduction of strong and fine-tuned
interactions and relies on the adjustments of parameters that can be
controlled in hybrid state-of-the-art Josephson junction platforms.
Additionally, we discussed numerical and experimental observables for the
detection of the main features of our system.

Our study poses the basis for further explorations of the correspondence
between nontrivial interacting conformal field theories and hybrid JJ systems
characterized by high harmonics terms. The ladder we devised, in particular,
provides a useful tool to engineer systems with exotic topological order in
two-dimensional systems: an array of these tricritical systems opens the way
to realize Fibonacci topological superconductors [20, 19] with universal non-
Abelian anyons.

#### Acknowledgements.-

We thank L. Banszerus, A. Cappelli, C. Marcus, G. Mussardo, C. Schrade and S.
Vaitiekėnas for fruitful discussions. We acknowledge support from the Deutsche
Forschungsgemeinschaft (DFG) project Grant No. 277101999 within the CRC
network TR 183 (subprojects B01 and C01). L.M. and M.B. are supported by the
Villum Foundation (Research Grant No. 25310). N.T. and M.R. are further
supported by the DFG under Germany’s Excellence Strategy – Cluster of
Excellence Matter and Light for Quantum Computing (ML4Q) EXC 2004/1 –
390534769. The authors gratefully acknowledge the Gauss Centre for
Supercomputing e.V. (www.gauss-centre.eu) for funding this project by
providing computing time through the John von Neumann Institute for Computing
(NIC) on the GCS Supercomputer JUWELS at the Jülich Supercomputing Centre
(JSC) (Grant NeTeNeSyQuMa) and the FZ Jülich for JURECA (institute project
PGI-8)[71].

## Supplemental materials for “Quantum simulation of the tricritical Ising
model in tunable Josephson junction ladders”

##  Appendix A Triple Josephson junction element

###  A.1 Higher harmonics expansion

In this section, we briefly analyze the decomposition of the energy-phase
relation of the triple JJ into harmonic terms $\mu_{n}$ that we introduced in
Eq. (2) of the main text. Assuming that each semiconducting/superconducting
junction is described by a single quantum channel, the potential of triple JJ
element

$V_{J}\left(\varphi\right)=-\Delta\left(\sqrt{1-T_{1}\sin^{2}\left(\frac{\varphi-\Phi_{1}}{2}\right)}+\sqrt{1-T_{2}\sin^{2}\left(\frac{\varphi\vphantom{\Phi}}{2}\right)}+\sqrt{1-T_{3}\sin^{2}\left(\frac{\varphi+\Phi_{2}}{2}\right)}\right),$
(S7)

can be expanded as $V_{J}=\sum_{n}\mu_{n}\cos{\left(n\varphi\right)}$, where
$\varphi$ is the SC phase difference. To maintain the reflection symmetry
$\varphi\to-\varphi$ we have imposed $\Phi_{1}=\Phi_{2}=\Phi$ and
$T_{1}=T_{3}$. The full expression of $\mu_{n}$ involves the elliptic
integrals

$\mu_{n}=\int_{-\pi}^{\pi}\dfrac{d\varphi}{\pi}\,V_{J}\left(\varphi\right)\cos{\left(n\varphi\right)},$
(S8)

whose analytical solution cannot be found. However, for small transparencies
$T_{i}\ll 1$, we can approximate them as follows:

$\begin{split}\mu_{1}/\Delta&=-\dfrac{1}{512}\left(T_{2}\left(128+32T_{2}+15T_{2}^{2}\right)+2T_{1}\left(128+32T_{1}+15T_{1}^{2}\right)\cos{\Phi}\right)+O\left(T_{i}^{4}\right)\\\
\mu_{2}/\Delta&=\dfrac{1}{256}\left(T_{2}^{2}\left(4+3T_{2}\right)+2T_{1}^{2}\left(4+3T_{1}\right)\cos{2\Phi}\right)+O\left(T_{i}^{4}\right)\\\
\mu_{3}/\Delta&=-\dfrac{1}{512}\left(T_{2}^{3}+2T_{1}^{3}\cos{3\Phi}\right)+O\left(T_{i}^{4}\right)\\\
\mu_{4}/\Delta&=O\left(T_{i}^{4}\right).\end{split}$ (S9)

In this limit, it is evident that the potential $V_{J}$ is mostly determined
by the first harmonic term $\cos{\varphi}$ with $\mu_{1}<0$, as long as the
flux $\Phi$ is such that $\cos{\Phi}>0$. Numerical evaluation of the integrals
(S8) shows that this is true also in the large transparencies limit.

The situation is different if we consider fluxes such that $\cos{\Phi}<0$. In
particular, one can fine-tune the external parameters to make $\mu_{1}$
vanish. Moreover, for $\Phi=2\pi/3$ and $T_{1}=T_{2}$ both $\mu_{1}$ and
$\mu_{2}$ vanish as a consequence of destructive interference of tunneling
events of one and two Cooper pairs through the three junctions. In this case
only triplet of Cooper pairs can jump between the two SC islands with
amplitude $|\mu_{3}|$. One can also check that, in the considered geometry,
the contribution $\mu_{4}$ is always at least one order of magnitude smaller
than the other terms as showed in Fig. S4. Therefore, given the ability of
controlling both the transparencies of the hybrid junctions through external
gates and the magnetic flux $\Phi$ piercing the two loops, we can tune
independently the ratios between the first three harmonics amplitudes in (S7).
In particular, the results discussed in the main text require that only the
transparencies of the external junctions, $T_{1}$ and $T_{3}$, need to be
tuned, whereas $T_{2}$ does not qualitatively affect the appearance of the
tricritical Ising point. This constitutes an advantage for experimental
realizations since the external junctions can more easily be controlled via
electrostatic gates.

Figure S4: The amplitudes of the first four harmonics $\mu_{n}$ with
$n=0,1,\ldots,4$ as a function of the triple JJ parameters. The white lines
mark the boundary of the symmetry-broken regime. We set $T_{2}=0.6$ and the SC
gap $\Delta$ induced by proximity in the semiconductors is fixed at
$\Delta=50$ in units of $E_{J}$ and does not influence the ratio between the
$\mu$ coefficients (S8).

###  A.2 Multichannel case

In the case of several transport channels in each of the junctions, the
Josephson energy-phase relation is given by the sum of the related
contributions:

$\mathcal{E}_{J}^{(p)}=-\sum_{i=1}^{M_{p}}\Delta\sqrt{1-T_{p}^{(i)}\sin^{2}{\left(\phi/2\right)}},$
(S10)

where $T_{p}^{(i)}$ represents the transparency of the $i$th channel in the JJ
$p$, and $M_{p}$ is the number of channels in the junction. For disordered
multichannel junctions, these transport coefficients $T_{p}^{(i)}$ follow a
bimodal distribution [72], with a few high-transparency channels resulting in
a nonsinusoidal current response. A complete generalization of our results to
the multichannel case goes beyond the scope of this supplemental section.
However, a qualitative analysis of its effects is needed. In particular, one
essential feature of our triple JJs element is the symmetry between the two
external junctions.

Experimental results in gate-tunable device showed that the nonsinusoidal
effects are overall well-approximated by one JJ with $M^{*}$ high-transparency
channels with the same average $T^{*}$, such that the current phase relation
reads [73, 11]

$I\left(\varphi\right)=\dfrac{e\Delta
M^{*}T^{*}}{\hbar}\dfrac{\sin{\left(\varphi\right)}}{\sqrt{1-T^{*}\sin^{2}{\left(\varphi/2\right)}}}.$
(S11)

Therefore, the non linear function in Eq. (1) in the main text well
approximate the energy-phase relation also in the multichannel case.

Moreover, in this approximation, one can assume that the external voltage gate
$V_{G}$ affects only the number of channels $M^{*}$ and not the average
transparency $T^{*}$ [11]. In this case, the symmetry between the external JJs
is lifted by a finite difference between the two average transparencies
$T^{*}_{1}-T^{*}_{3}\neq 0$, which is almost independent of the voltage gates
$V_{G1}$ and $V_{G3}$. However, tuning the number of open channels $M^{*}_{1}$
and $M^{*}_{3}$ via the voltage gates provides a way to mitigate this explicit
symmetry breaking. Finally, potential asymmetries in the magnetic fluxes cause
a splitting in energy of the minima of the potential $V_{J}$ which is linear
in $\Phi_{1}-\Phi_{3}$. However, this effect can also be used to mitigate the
asymmetry caused by the mismatch of the transparencies $T_{1}\neq T_{3}$ and
restore the degeneracy of the minima of $V_{J}$.

##  Appendix B Ladder: further details

###  B.1 Staggered magnetic fluxes

Interacting bosons on a ladder with uniform magnetic fields exhibit are
characterized by the onset of several chiral many-body phases, including the
Meissner phase. For our purposes the onset of the Meissner effect may be
detrimental, because it breaks the emergent Lorentz invariance in the QFT and
may compete with the phases and critical points discussed in the main text.

Additionally, to obtain a quantum simulation of the three-frequency sine-
Gordon model, each rung triple JJ must be characterizes by the same $V_{J}$.
This condition is, in the general case, fulfilled only by staggered patterns
of magnetic fluxes.

We present two viable flux configurations which are schematically represented
in Fig. S5(a) and (b). The solution (a) relies on the parity property of the
local potential $V_{J}$ under $\Phi\to-\Phi$ and enables the engineering of a
ladder geometry where the magnetic flux within two subsequent rungs $\Phi_{\rm
int}$ vanishes. This preserves time-reversal invariance in the effective QFT.
However, this approach leads to the experimental challenge of controlling
nonuniform magnetic fields along the ladder.

A convenient construction to realize the configuration (a) in experimental
devices is depicted in Fig. S5(c). To stagger the magnetic fluxes within two
subsequent triple JJ elements, we design the ladder in a ’snake’ configuration
and control the magnetic field by introducing a current $I_{\rm ext}$ through
the line schematically represented in Fig. S5. Alternatively, a local control
of multiple fluxes can be achieved with the techniques adopted by modern
quantum processors based on transmon qubits [74].

An alternative flux configuration, Fig. S5(b) results in the same potentials
$V_{J}$ on each rung and relies on compensating the magnetic fluxes of the
triple JJs with opposite fluxes in the ladder plaquettes, thus setting
$\Phi_{\rm int}=-2\Phi$ between each rung. The possibility of introducing
additional integer fluxes in each loop, thus replacing $\Phi_{\rm
int}\to\Phi_{\rm int}+2\pi$ may also offer an alternative to implement the
configuration (b) with uniform magnetic fluxes. To tune the system at the
tricritical point in this scenario, however, it is required to known a priori
the parameter $T_{2}$ of the ladder, such that the correct magnetic fluxes can
be set by designing superconducting circuits of a suitable area.

Figure S5: (a) and (b) illustrate the two configurations of nonuniform
magnetic fluxes that allow to engineer the potential $V_{J}$ on the ladder and
fulfill the required $\mathbb{Z}_{2}$ symmetry. In configuration (a), the
fluxes are staggered within consecutive triple Josephson junction (JJ)
elements, while in (b), a $\Phi_{\rm int}=-2\Phi$ flux is inserted in the
plaquettes of the effective ladder. Panel (c) shows the physical realization
of configuration (a) achieved through a snake geometry and the insertion of a
line with a tunable external current $I_{\rm ext}$.

###  B.2 Bosonization

In this section, we will review the main steps of the connection between the
lattice Hamiltonian in (4) in the main text and the three-frequency sine-
Gordon quantum field theory. At low temperature $K_{B}T<\Delta_{c}$ each SC
island of our lattice corresponds to a condensate of $N_{c}$ Cooper pairs with
gap $\Delta_{c}$ and a well defined complex order parameter, the SC phase
$\hat{\varphi}_{\alpha,j}$. The residual charge around $N_{c}$ is represented
by the operator $\widehat{N}_{\alpha,j}$ dual to the SC phase. In the long
wavelength limit, we can use an effective continuum description in terms of
the Bose fields $\hat{\theta}_{\alpha}(x)$ and $\hat{\varphi}_{\alpha}(x)$
[33], fulfilling commutation relations:

$\left[\hat{\theta}_{\alpha}(y),\hat{\varphi}_{\beta}(x)\right]=-i\pi\delta_{\alpha\beta}\Theta\left(y-x\right)\,,$
(S12)

where $\Theta$ indicates the Heaviside step function. The weak interactions
case $E_{C},\;V_{\perp},\;\ll E_{J}$ we considered allows us to neglect fast-
oscillating contributions in the Cooper-pair density and write
$N_{\alpha,j}\approx-a\dfrac{\partial_{x}\theta_{\alpha}(x)}{\pi}$, with
$j=xa$. In the harmonic approximation for the Josephson interaction along the
legs, the low-energy lattice Hamiltonian can be written as

$\hat{H}=\sum_{\alpha={a},{b}}\left[-E_{J}\int
dx\;a\left(\partial_{x}\hat{\varphi}_{\alpha}\left(x\right)\right)^{2}+\dfrac{E_{C}a}{\pi^{2}}\int
dx\;\left(\partial_{x}\hat{\theta}_{\alpha}\left(x\right)\right)^{2}\right]+\sum_{n=1}^{3}\frac{\mu_{n}}{a}\int
dx\;\cos{\left(n\left(\hat{\varphi}_{{a}}-\hat{\varphi}_{{b}}\right)\right)}.$
(S13)

By rotating the fields
$\hat{\varphi}_{c/s}(x)=\left(\hat{\varphi}_{{a}}(x)\pm\hat{\varphi}_{{b}}(x)\right)/\sqrt{2}$
and the corresponding dual ones $\hat{\theta}_{c/s}(x)$, we obtain the
Hamiltonian (5) in the main text with the perturbative relations

$K_{c/s}=\pi\sqrt{\dfrac{E_{J}}{\left(2E_{c}\pm
V_{\perp}\right)}}\qquad\text{and}\qquad u_{c/s}=a\sqrt{E_{J}\left(2E_{C}\pm
V_{\perp}\right)}.$ (S14)

In general, a finite intra-leg capacitance $C_{L}$ among adjacent islands
leads to a long range interaction stemming from the inverse capacitance matrix
[52] with screening length $\lambda=a\sqrt{C_{L}/C_{g}}$, where $C_{g}$ is the
self capacitance. However, this may be ignored as long as one is interested in
the physics of modes with energies lower than $u_{c/s}/\lambda$.

From a perturbative point of view the plasma frequency of the spin sector
$u_{s}/a=\Lambda\simeq\sqrt{E_{J}\left(2E_{c}-V_{\perp}\right)}$ defines a UV
cut-off that allows us to define the dimensionless coupling
$\tilde{\mu}_{n}=\mu_{n}/\Lambda$ in the sine-Gordon Euclidean action,

$S\left[\varphi_{s}(x,\tau)\right]=\dfrac{1}{2\pi}\int
dxd\tau\;K_{s}\left(\left(\partial_{\tau}\varphi_{s}\right)^{2}+\left(\partial_{x}\varphi_{s}\right)^{2}\right)-\sum_{n=1}^{3}\frac{\tilde{\mu}_{n}}{a^{2}}\int
dxd\tau\;\cos{\left(\sqrt{2}n\varphi_{s}\right)},$ (S15)

where we have rescaled the imaginary time $\tau\to u_{s}\tau$. The operators
$\widehat{\mathcal{O}}_{n}=\cos{\left(\sqrt{2}n\hat{\varphi}_{s}\right)}$
correspond to primaries of the unperturbed free boson $c=1$ theory with
scaling dimensions

$\Delta_{n}=\dfrac{n^{2}}{2K_{s}},$ (S16)

and they drive the LL to a massive phase, namely they are relevant, only when
$\Delta_{n}<2$.

Note that the charge sector remains massless as there is no sine-Gordon
potential for $\hat{\varphi}_{c}$. We checked the validity of this statement
in our lattice simulation. In the LL liquid phase the density correlation
functions is expected to show the following power-law decay

$\braket{\widehat{\rho}_{\rm tot}(x)\widehat{\rho}_{\rm
tot}(y)}\sim\frac{2}{\pi^{2}}\left\langle\partial_{x}\theta_{c}(x,\tau)\;\partial_{y}\theta_{c}(y,\tau)\right\rangle=\frac{K_{c}}{\pi^{2}}\frac{1}{\left|x-y\right|^{2}}.$
(S17)

In the ladder model, the operator $\widehat{\rho}_{\rm tot}(x)$ corresponds to
the total rung density offset $\widehat{N}_{{\rm
tot},j}-\braket{\widehat{N}_{{\rm tot},j}}$ with $\widehat{N}_{\rm
tot}=\widehat{N}_{{a},j}+\widehat{N}_{{b},j}$. We explicitly checked the decay
of Eq. (S17) for each point of the phase diagram by fitting a power-law decay
[Fig. S6]. The so found $K_{c}$ parameters are in a good agreement with the
perturbative approximations given by Eq. (S14). This confirms the validity of
the field theoretical approach in the low energy regime of the ladder.

Figure S6: Connected part of the two-point correlation function of the total
density operator
$\widehat{N}_{\text{tot},j}=\left(\widehat{N}_{{a},j}+\widehat{N}_{{b},j}\right)$
taken at a random position in the phase-diagram at $h_{1}\approx-0.3$ and
$h_{2}\approx 0.47$. The red line is the result of a fit by a function
$f(j)=K_{c}/\pi^{2}j^{-\alpha}$. The exponent $\alpha=2$ matches the
predictions of bosonization; the obtained Luttinger parameter is close to the
predicted value: $K_{c}^{\rm pert}\approx 2.61$.

On the other hand, the spin sector (S15) is subject to the different relevant
interactions in Eq. (S15) which tend to order the SC phase difference
$\hat{\varphi}_{s}$. In Ref. [44] the author shows that this quantum field
theory flows to a tricritical Ising point with central charge $c=7/10$ for
suitable values of the coupling constants $\mu$. Despite the absence of any
non-perturbative mappings between our lattice operators and the massless
excitations of this field theory, we can exploit the Ginzburg-Landau
representation of the TCI CFT to gain insight about this relation.

The operator content of the CFT is split in the odd and even sector with
respect to the $\mathds{Z}_{2}$-symmetry and is characterized by 6 primary
fields: the identity $I$, 4 relevant operators
$\sigma,\;\epsilon\;\sigma^{\prime},\;\epsilon^{\prime}$ ($\Delta<2$) and one
irrelevant $(\Delta>2)$ operator. The Ginzburg-Landau Lagrangian
representation of the TCI corresponds to [75]

$\mathcal{L}=\dfrac{K_{s}}{2\pi}\varphi_{s}\left(\partial_{x}^{2}+\dfrac{\partial_{\tau}^{2}}{u_{s}^{2}}\right)\varphi_{s}-\lambda_{2}:\varphi_{s}^{2}:-\lambda_{4}:\varphi_{s}^{4}:-\lambda_{6}:\varphi_{s}^{6}:,$
(S18)

where $::$ indicates the normal ordering with respect to the tricritical point
CFT. In the mean-field limit $K_{s}\gg 1$, we can build an approximate mapping
bewteen local operators in our theory and the primary fields (see also Ref.
[27]),

$\displaystyle\varphi_{s}(x)\to\sigma(x),\quad\left(h_{\sigma},\bar{h}_{\sigma}\right)=\left(\dfrac{3}{80},\dfrac{3}{80}\right)$
(S19)
$\displaystyle:\varphi_{s}^{2}(x):\to\epsilon(x),\quad\left(h_{\epsilon},\bar{h}_{\epsilon}\right)=\left(\dfrac{1}{10},\dfrac{1}{10}\right)$
$\displaystyle:\varphi_{s}^{3}(x):\to\sigma^{\prime}(x),\quad\left(h_{\sigma^{\prime}},\bar{h}_{\sigma^{\prime}}\right)=\left(\dfrac{7}{16},\dfrac{7}{16}\right)$
$\displaystyle:\varphi_{s}^{4}(x):\to\epsilon^{\prime}(x),\quad\left(h_{\epsilon^{\prime}},\bar{h}_{\epsilon^{\prime}}\right)=\left(\dfrac{3}{5},\dfrac{3}{5}\right),$

which implies the expansion of the local order operator $\hat{J}_{\perp}$ in
terms of the most relevant operator $\sigma$ close to the critical point,

$\hat{J}_{\perp}(x)=\sin{\left(\sqrt{2}\hat{\varphi}_{s}(x)\right)\sim\hat{\varphi}_{s}(x)+\ldots\to\sigma(x)+\ldots}$
(S20)

In the previous expansion the dots indicate less relevant operator
contributions.

##  Appendix C Charge basis

For the numerical simulations, we formulated the Hamiltonian (4) from the main
text in the charge basis. In this basis the operator $\widehat{N}_{\alpha,j}$
is diagonal and defines how the number of Cooper pairs differs from the
average occupation on the island $(\alpha,j)$:

$\widehat{N}_{\alpha,j}=\mathrm{diag}\left(\dots,-2,-1,0,1,2,\dots\right)\,.$
(S21)

Using this choice, it is easy to show that $e^{i\hat{\varphi}_{\alpha,j}}$
must to be of the form

$e^{i\hat{\varphi}_{\alpha,j}}=\begin{pmatrix}&\ddots&&&&&&\\\ &&0&1&&&&\\\
&&&0&1&&&\\\ &&&&0&1&&\\\ &&&&&\ddots&&\\\
\end{pmatrix}_{\alpha,j}\equiv\widehat{\Sigma}_{\alpha,j}^{-}$ (S22)

for the commutator $[\widehat{N},\widehat{\Sigma}^{-}]=-\widehat{\Sigma}^{-}$
to hold. Further, in order to represent these operators in our simulations, we
have to truncate the number of possible charge states

$\widehat{N}_{\alpha,j}=\mathrm{diag}\left(-N_{\rm max}\dots,-2,-1,0,1,2,\dots
N_{\rm max}\right)\,,$ (S23)

i.e. we adopts a truncated local Hilbert-space of dimension $2N_{\rm max}+1$
per each SC island. We can control the error caused by this truncation by
varying $N_{\rm max}$ until we reach convergence in all observables.
Alternatively, we can measure the probability
$\braket{\hat{P}^{n}_{\alpha,j}}$ of finding an excitation $n$ on the island
$(\alpha,j)$. By ensuring that $N_{\rm max}$ is large enough to have
negligible weight $\braket{\hat{P}^{N_{\rm max}}_{\alpha,j}}<\epsilon$ we can
claim to be converged in $N_{\rm max}$. In practice we found that $N_{\rm
max}=8$ gives $\braket{\hat{P}^{N_{\rm max}}_{\alpha,j}}\sim 10^{-9}$. The
Hamiltonian used for the simulation finally reads
$\widehat{H}=\sum\limits_{j=0}^{L}\widehat{h}_{j,j+1}$ with:

$\begin{split}\hat{h}_{j,j+1}=&\sum_{\alpha=a,b}\left[E_{c}\left(\widehat{N}_{\alpha,j}\right)^{2}-\frac{E_{J}}{2}\left(\widehat{\Sigma}_{\alpha,j}^{+}\widehat{\Sigma}_{\alpha,j+1}^{-}+\widehat{\Sigma}_{\alpha,j}^{-}\widehat{\Sigma}_{\alpha,j+1}^{+}\right)\right]\\\
{}&+V\widehat{N}_{a,j}\widehat{N}_{b,j}+\frac{\mu_{1}}{2}\left(\widehat{\Sigma}_{a,j}^{+}\widehat{\Sigma}_{b,j}^{-}+\widehat{\Sigma}_{b,j}^{+}\widehat{\Sigma}_{a,j}^{-}\right)\\\
{}&+\frac{\mu_{2}}{2}\left(\left(\widehat{\Sigma}_{a,j}^{+}\right)^{2}\left(\widehat{\Sigma}_{b,j}^{-}\right)^{2}+\left(\widehat{\Sigma}_{b,j}^{+}\right)^{2}\left(\widehat{\Sigma}_{a,j}^{-}\right)^{2}\right)\\\
{}&+\frac{\mu_{3}}{2}\left(\left(\widehat{\Sigma}_{a,j}^{+}\right)^{3}\left(\widehat{\Sigma}_{b,j}^{-}\right)^{3}+\left(\widehat{\Sigma}_{b,j}^{+}\right)^{3}\left(\widehat{\Sigma}_{a,j}^{-}\right)^{3}\right)\end{split}$
(S24)

##  Appendix D Further numerical results

In this section, we present additional numerical indications about the
different nature of the transitions across the phase diagram.

###  D.1 Hysteresis and gap jump at the first-order transition

First of all, we present additional evidence of first-order phase transitions
(FOPTs) along the horizontal cuts at $X_{2}=0.52$ (between the disordered
phase I and the ordered phase II) and at $X_{2}=0.6$ (between phases I and
III).

One significant indicator involves the distinct behavior of the lowest energy
excitation in the spin sector. Its energy corresponds to the system’s gap,
which can be extracted (see section below) from the transfer matrix spectrum
as shown in Fig. S7. By following the corresponding eigenvalue of the transfer
matrix $\lambda_{1}$, we can extract the gap of the spin sector
$\Delta_{s}=-\log{\lambda_{1}}$. Across a second-order phase transition, the
physical gap closes and, in the numerical VUMPS simulations, this is marked by
a minimum in $\Delta_{s}$ [panel (c)] which approaches zero by increasing the
bond dimension. Across a FOPT, instead, the spin gap remains finite [panels
(a) and (b)]. Panels (a) and (b) respectively depict the typical behaviors of
the FOPT between the two disordered phases and between phase II and phase I.
In the latter case, the related order parameter displays a very weak
variation, resulting in an almost continuous behavior of $\Delta_{s}$.

This behavior is reflected also in the analysis of the hysteresis in the order
parameter and the many-body ground state energy, as illustrated in Fig. S8.

A discontinuity in the first derivative of the energy density is observed in
the FOPT cases, which is absent in the second-order transition at $X_{2}=0$
and indicates the crossing of the lowest energy levels. Furthermore, by
altering the minimization procedure at each point $X_{1}$ and initializing the
ground state with the result from $X_{1}\pm\delta$, the variational algorithm
follows the corresponding branch, even within the opposite phase. This can be
interpreted as a hysteresis effect induced by the orthogonality of these two
states around the crossing point.

Also in this case the features of the FOPT are stronger between the two
disordered phases – panel S8 (b) is depicted with a magnified energy scales
with respect to panel (a). The discontinuity of the derivative
$\partial\varepsilon/\partial X_{1}$ is around $29.5\;E_{J}$ in panel (a) and
$21.7\;E_{J}$ in panel (b). This is physically related to the jump of the
average loop current circulating around each triple JJs element, namely
$\hat{J}_{\rm loop}=\partial\hat{H}/\partial\Phi$.

Figure S7: The gap in the spin sector $\Delta_{s}$ is determined by tracking
the second largest eigenvalue, $\lambda_{1}$, of the transfer matrix within
the spin sector. The results are obtained for the two cuts shown in Fig. 3a)
and an additional line between phases I and III: a) $X_{2}=0.6$, b)
$X_{2}=0.52$ and c) $X_{2}=0$. The red points correspond to the gap
$\Delta_{s}=-\log(\lambda_{1})$ which remains finite across the FOPTs in panel
(a) and (b), while displaying the gap-closing feature of a second-order phase
transition in panel (c). Figure S8:  Energy density $\varepsilon$ of the
groundstate obtained along the same lines of Fig. S7: a) $X_{2}=0.6$, b)
$X_{2}=0.52$ and c) $X_{2}=0$. The red triangles in the case of a) and b) are
obtained by minimizing the Hamiltonian $H({\bf X}+\delta)$ by starting from
one of the two groundstates left/right of the meeting point of the two
branches. The minimization procedure follows these branches instead of falling
into the true ground state. The absence of such an effect for $X_{2}=0$ is
another indication for a FOPT in the case a) and b), but a second order phase
transition for c).

###  D.2 Scaling and critical exponents Ising phase transition

Figure S9:  Extrapolation of the critical exponent $\nu$ from the scaling
features of the fidelity at $X_{2}=0$. (a): Two-parameter fit of the relation
$\mathcal{F}\left(X_{1c},\delta\right)=a\left|\delta\right|^{\nu}$ at bond
dimension $D=1000$. The extracted values of $\nu$ increase with the increasing
bond dimension (inset). (b): Fit of the fidelity susceptibility
$\chi_{\mathcal{F}}=b\left|X_{1}-X_{1c}\right|^{\nu-2}$ away from the critical
point with a fixed $\delta\ll\left|X_{1}-X_{1c}\right|$. The plot in log-log
scale is shown in the inset. The position of the critical point $X_{1c}$ is
obtained from the collapse of the spin correlation length.

In this subsection, we focus on characterizing the critical exponents $\nu$
and $\beta$, which describe how the correlation length diverges and the order
parameter approaches zero across the continuous phase transition at the
$X_{2}=T_{1}\sin(\Phi)=0$ cut. As highlighted in the main text, the measured
values indicate that the transition belongs to the Ising universality class
with $\nu_{\rm IS}=1$ and $\beta_{\rm IS}=1/8$. To extract these exponents, we
relied on scaling properties of three different quantities: the log-fidelity
per site $\mathcal{F}$ (and its susceptibility $\chi_{\mathcal{F}}$), the
correlation length of the spin sector $\xi_{s}$ and the order parameter
$\hat{J}^{(2e)}_{\perp}$.

We determine the critical exponent $\nu$ through two different methods based
on the fidelity scaling, both yielding values near $\nu_{\rm IS}=1$ [Fig. S9].
The first approach involves fitting the non-analytic behavior of the log-
fidelity per site at the critical point, showing a consistent increase towards
$\nu=1$ as the bond dimension $D$ grows [Fig S9 inset]. Alternatively, by
analyzing the divergence pattern of the fidelity susceptibility along a
horizontal cut, we obtain $\nu=1.00(3)$.

To take into account finite bond dimension corrections, we employed the finite
entanglement scaling discussed in [58] for the spin correlation length
$\xi_{s}$. Similarly to finite size effects, the finite bond dimension
introduces an artificial length scale rendering all correlation functions
exponential decaying even at critical points. This can be understood by a
relevant perturbation of the underlying CFT. However, in the $D\to\infty$
limit, the gapless nature of the model must be restored. This artificial
length scale is associated with the critical exponent $\kappa$:

$\xi_{D}\sim D^{\kappa}$

and we use this relation to define the following scaling ansatz [58]

$\xi_{D}=D^{\kappa}f\left(D^{\frac{\kappa}{\nu}}\frac{|X_{1}-X_{1c}|}{X_{1c}}\right)\,,\quad
f(x)\sim\begin{cases}\text{const}\ ,\,&x\to 0\\\ \frac{1}{x^{\nu}}\ ,\,&x\gg
1\end{cases}$ (S25)

where $\nu$ is the critical exponent of the correlation length in the infinite
bond dimension case. We use this ansatz to determine the critical point
$X_{1c}$ and to extract the critical exponent $\nu$ discussed in the main
text.

Additionally, to extract the critical exponent $\beta$ we employ the scaling
of the expectation value of the single-particle current
$\hat{J}^{(2e)}_{\perp}$ close to the critical point. Indeed, this operator
plays the role of the Ising magnetization which is odd under the
$\mathds{Z}_{2}$-symmetry $\hat{\varphi}_{s}\to-\hat{\varphi}_{s}$. By fitting
the expected scaling behaviour $\left|X_{1}-X_{1c}\right|^{\beta}$, we obtain
the critical exponent $\beta=0.125(3)$ [Fig. S10].

These results collectively indicate that our findings concerning the
transition from the ordered to the disordered phase sufficiently far from the
first order discontinuities are compatible with the Ising universality class
with $\nu_{\rm IS}=1$ and $\beta_{\rm IS}=1/8$.

The critical exponents $\kappa$ extracted for the spin correlation length at
the second order transitions are typically smaller than one. This implies a
huge increase of the bond dimension in order to faithfully capture the
algebraic decay of correlation functions over a long distance. Taking the
example of the $X_{2}=0$ cut from the main text with $\kappa\approx 0.8$. The
largest correlation length obtained for $X_{2}$ is $\xi_{s}\approx 30$ for a
bond dimension of $D=1000$. Using the scaling behavior $\xi_{s}\sim D^{0.8}$
we estimate that a bond dimension $D^{\star}\approx 4500$ is necessary to get
$\xi_{s}\approx 100$ sites, and $D^{\star}\approx 18000$ for $\xi_{s}\approx
300$ sites.

Figure S10: The single-particle current $J^{(2e)}_{\perp}$ plays the role of
the effective magnetization at the Ising critical point, displaying a scaling
behavior $|X_{1}-X_{1c}|^{\beta}$ with the fitted value $\beta=0.125\pm 0.003$
(red curve). The critical point $X_{1c}$ is fixed by the collapse of $\xi_{s}$
obtained by using Eq. (S25). The discrepancy with the numerical points is due
to finite entanglement effects that shifts the position of the critical point
at finite bond dimensions.

###  D.3 Central charge

Figure S11: The color plot shows the distance from 1 of the total central
charge $c$ obtained from the entanglement entropy. Away from the transition
between the phases I and II, we detect the massless charge sector with $c=1$
(blue areas). On the other hand, along the second order transition line, the
total central charge fluctuates between $c\sim 1.5-1.8$ until it starts
decreasing to 1 where the second-order line is replaced by the the FOPT.

Given the separation of the two sectors in our model, the entanglement entropy
of the system is predicted to display a typical divergence
$S=c_{c}/6\log(\xi_{c})+c_{s}/6\log(\xi_{s})$ [76] in proximity of the second-
order phase transition. However, at the phase transition the ratio between the
two correlation length remains constant when increasing the bond dimension.
Therefore we can consider the simpler expression $S=c/6\log(\xi)$ where $\xi$
is a diverging correlation length and $c$ represents the central charge of the
complete underlying CFT.

Therefore, to investigate the appearance of the TCI point, we focus on the
zoomed region shown in Fig S11 and determine the total central charge $c$ of
the system by fitting this logarithmic divergence. The correlation length
$\xi_{D}$ of the critical system is obtained from the second largest
eigenvalue of the transfer matrix evaluated for different values of the bond
dimension.

The strong finite entanglement effects have a quantitative impact on our
results. In particular, they lead to strong fluctuations that result in noisy
estimates that do not allow us to clearly observe the predicted non-monotonic
behavior from $c=1+1/2$ to $c=1+7/70$ [Fig. S11]. Nevertheless, following the
transition, we detect an expected drop of the central charge in correspondence
of the FOPT line, which confirms the connection between first and second-order
phase transitions.

##  Appendix E Extraction Of Correlation Lenghts

All numerical results presented in this latter are obtained by the VUMPS
algorithm presented in Ref. [36]. The concrete implementation uses the ITensor
library [77]. This ansatz operates directly in the thermodynamic limit by
enforcing translational invarance. The class of ansatz states is characterized
by the set of matrices $\\{A_{L}^{\sigma},A_{C}^{\sigma},A_{R}^{\sigma}\\}$,
with $\sigma$ enumerating the physical local states. From this set of
matrices, the state $\ket{\psi}$ is represented as

$\ket{\psi}=\sum_{\\{\sigma\\}}\operatorname{Tr}\left[\dots
A_{L}^{\sigma_{j-2}}A_{L}^{\sigma_{j-1}}A_{C}^{\sigma_{j}}A_{R}^{\sigma_{j+1}}A_{L}^{\sigma_{j+2}}\dots\right]\ket{\dots\sigma_{j-2}\sigma_{j-1}\sigma_{j}\sigma_{j+1}\sigma_{j+2}\dots}\,.$

The matrices $A_{L}^{\sigma}$ and $A_{R}^{\sigma}$ fulfill
$\sum_{\sigma}(A_{L}^{\sigma})^{\dagger}A_{L}^{\sigma}=\sum_{\sigma}A_{R}^{\sigma}(A_{R}^{\sigma})^{\dagger}=\operatorname{\mathds{1}}$
and special equivariance relations to ensure the translational invariance of
the ansatz, see Fig. S12.

$\ket{\psi}=$$\dots$$A_{\rm L}$$A_{\rm C}$$A_{\rm
R}$$\sigma_{j-1}$$\sigma_{j}$$\sigma_{j+1}$$\dots$$A_{\rm L}$$A_{\rm
C}$$=$$A_{\rm C}$$A_{\rm R}$$a)$$b)$ Figure S12: (a): VUMPS ansatz in the
central gauge. (b): Equivariance property to ensure translational invariance

Using the transfer-matrix of the system, defined by

$\mathcal{T}_{L}\coloneqq\sum_{\sigma}A_{L}^{\sigma}\otimes\bar{A}_{L}^{\sigma}\,,$
(S26)

and the two transfer-matrices with operator insertion

$\mathcal{T}_{L}^{O}\coloneqq\sum_{\sigma,\tau}O_{\sigma,\tau}A^{\sigma}_{L}\otimes\bar{A}^{\tau}_{L}\,,\quad\mathcal{T}_{C}^{K}\coloneqq\sum_{\sigma,\tau}K_{\sigma,\tau}A^{\sigma}_{C}\otimes\bar{A}^{\tau}_{C}\,,\quad$
(S27)

where $\bar{z}$ denotes the complex conjugation of $z$, one can represent the
correlation function of two arbitrary operators $\hat{O}$ and $\hat{K}$ as,
Fig. S13:

$\begin{split}\braket{\hat{O}_{j}\hat{K}_{j+l}}&=\bra{\operatorname{\mathds{1}}}\mathcal{T}^{O}_{\rm
L}\left(\mathcal{T}_{\rm L}\right)^{l-1}\mathcal{T}^{K}_{\rm
C}\ket{\operatorname{\mathds{1}}}=\sum_{n\geq
0}\lambda_{n}^{l-1}\alpha^{O}_{n}\,\beta^{K}_{n}=\sum_{n\geq
0}e^{-\frac{l-1}{\xi_{n}}}c^{O,K}_{n}\\\
\alpha^{O}_{n}&=\braket{\operatorname{\mathds{1}}}{\mathcal{T}_{O}}{R_{n}},\
\beta^{K}_{n}=\braket{L_{n}}{\mathcal{T}_{K}}{\operatorname{\mathds{1}}},\
\xi_{n}=-\frac{1}{\log(\lambda_{n})}\,.\end{split}$ (S28)

The second line in Eq. S28 is obtained after using the eigen decomposition of
the transfer-matrix

$\mathcal{T}_{L}=\sum_{n\geq
0}\lambda_{n}\ket{R_{n}}\bra{L_{n}}\,,\quad\braket{L_{n}}{R_{m}}=\delta_{m,n}\,.$
(S29)

$\braket{\hat{O}_{j}\hat{K}_{j+l}}\coloneqq$$\dots$$A_{\rm L}$$\bar{A}_{\rm
L}$$A_{\rm L}$$\bar{A}_{\rm L}$$A_{\rm L}$$\bar{A}_{\rm L}$$O$$\dots$$A_{\rm
L}$$\bar{A}_{\rm L}$$A_{\rm C}$$\bar{A}_{\rm C}$$K$$A_{\rm R}$$\bar{A}_{\rm
R}$$\dots$$l-1$$=$$A_{\rm L}$$O$$\bar{A}_{\rm L}$$A_{\rm C}$$\bar{A}_{\rm
C}$$K$ $\mathcal{T}_{\rm L}$ $l-1$ Figure S13: Correlation function in the
infinite system.

Using Eq. S28, it is straightforward to extract the asymptotic behavior of any
correlation function

$\braket{\hat{O}_{j}^{\phantom{{\dagger}}}\hat{K}_{j+l}^{\dagger}}\approx
c_{n^{\star}}^{O,K}\,e^{\frac{l}{\xi_{n^{\star}}}}+c_{0}^{O,K}\,.$

where $n^{\star}$ is the first $n>0$ in the descending sequence
$\lambda_{0}>|\lambda_{1}|\geq|\lambda_{2}|\dots$ with a non-zero operator
weight $c_{n}^{O,K}$ (assuming $\lambda_{n^{\star}}$ to be unique). The
contribution $c_{0}^{O,K}$ equals the product of expectation values
$\braket{\hat{O}_{j}}\braket{K^{\dagger}_{j}}$. In the case of
$\hat{O}=\hat{K}$ this asymptotic behavior can be used to extract the smallest
energy gap in the excitation spectrum generated by the operator $\hat{O}$. In
the main text, we applied this analysis to the current operator

$\hat{O}=\widehat{J}^{(2e)}_{\perp}\coloneqq\frac{i}{2}\left(\Sigma_{{a}}^{+}\Sigma_{{b}}^{-}-\Sigma_{{b}}^{+}\Sigma_{{a}}^{-}\right)\,.$

which can be interpreted as the magnetization order parameter in the field
theory $\sin{\left(\hat{\varphi}_{s}(x)\right)}$ odd under the
$\varphi_{s}(x)\to-\varphi_{s}(x)$ symmetry transformation. Thus,
$\hat{J}^{(2e)}_{\perp}$ is naturally associated to excitations in the spin-
sector exclusively. Very similar, one can extract the density of the
logarithmic fidelity $\mathcal{F}$ in the thermodynamic limit from the mixed
transfer-matrix

$\mathcal{T}_{L}^{\phi,\psi}\coloneqq\sum_{\sigma}A_{L}^{\phi,\sigma}\otimes\bar{A}_{L}^{\psi,\sigma}\,,$
(S30)

where $A_{L}^{\phi}$ defines the state $\ket{\phi}$ and $A_{L}^{\psi}$ the
state $\ket{\psi}$. Define $\lambda_{0}$ the smallest in maginutde eigenvalue
of $\mathcal{T}_{L}^{\phi,\psi}$, it is straigthforward to show:

$\mathcal{F}\coloneqq-\lim_{N\to\infty}\frac{1}{N}\log\left(\braket{\psi}{\phi}\right)=-\log(|\lambda_{0}|)\,.$

## References

  * Krogstrup  _et al._ [2015] P. Krogstrup, N. L. B. Ziino, W. Chang, S. M. Albrecht, M. H. Madsen, E. Johnson, J. Nygård, C. M. Marcus, and T. S. Jespersen, Epitaxy of semiconductor–superconductor nanowires, Nat. Mater. 14, 400 (2015).
  * Shabani  _et al._ [2016] J. Shabani, M. Kjaergaard, H. J. Suominen, Y. Kim, F. Nichele, K. Pakrouski, T. Stankevic, R. M. Lutchyn, P. Krogstrup, R. Feidenhans’l, S. Kraemer, C. Nayak, M. Troyer, C. M. Marcus, and C. J. Palmstrøm, Two-dimensional epitaxial superconductor-semiconductor heterostructures: A platform for topological superconducting networks, Phys. Rev. B 93, 155402 (2016).
  * Bøttcher  _et al._ [2018] C. G. L. Bøttcher, F. Nichele, M. Kjaergaard, H. J. Suominen, J. Shabani, C. J. Palmstrøm, and C. M. Marcus, Superconducting, insulating and anomalous metallic regimes in a gated two-dimensional semiconductor–superconductor array, Nat. Phys. 14, 1138 (2018).
  * Bøttcher  _et al._ [2022a] C. G. L. Bøttcher, F. Nichele, J. Shabani, C. J. Palmstrøm, and C. M. Marcus, The Berezinskii-Kosterlitz-Thouless Transition and Anomalous Metallic Phase in a Hybrid Josephson Junction Array, arXiv 10.48550/arXiv.2210.00318 (2022a), 2210.00318 .
  * Bøttcher  _et al._ [2022b] C. G. L. Bøttcher, F. Nichele, J. Shabani, C. J. Palmstrøm, and C. M. Marcus, Dynamical vortex transitions in a gate-tunable Josephson junction array, arXiv 10.48550/arXiv.2212.08651 (2022b), 2212.08651 .
  * Kuzmin  _et al._ [2019] R. Kuzmin, R. Mencia, N. Grabon, N. Mehta, Y.-H. Lin, and V. E. Manucharyan, Quantum electrodynamics of a superconductor–insulator phase transition, Nat. Phys. 15, 930 (2019).
  * Mukhopadhyay  _et al._ [2022] S. Mukhopadhyay, J. Senior, J. Saez-Mollejo, D. Puglia, M. Zemlicka, J. Fink, and A. P. Higginbotham, Superconductivity from a melted insulator, arXiv 10.48550/arXiv.2210.06508 (2022), 2210.06508 .
  * Cedergren  _et al._ [2017] K. Cedergren, R. Ackroyd, S. Kafanov, N. Vogt, A. Shnirman, and T. Duty, Insulating Josephson Junction Chains as Pinned Luttinger Liquids, Phys. Rev. Lett. 119, 167701 (2017).
  * Kjaergaard  _et al._ [2017] M. Kjaergaard, H. J. Suominen, M. P. Nowak, A. R. Akhmerov, J. Shabani, C. J. Palmstrøm, F. Nichele, and C. M. Marcus, Transparent semiconductor-superconductor interface and induced gap in an epitaxial heterostructure Josephson junction, Phys. Rev. Appl. 7, 034029 (2017).
  * Casparis  _et al._ [2018] L. Casparis, M. R. Connolly, M. Kjaergaard, N. J. Pearson, A. Kringhøj, T. W. Larsen, F. Kuemmeth, T. Wang, C. Thomas, S. Gronin, G. C. Gardner, M. J. Manfra, C. M. Marcus, and K. D. Petersson, Superconducting gatemon qubit based on a proximitized two-dimensional electron gas, Nature Nanotechnology 13, 915 (2018).
  * Ciaccia  _et al._ [2023] C. Ciaccia, R. Haller, A. C. C. Drachmann, T. Lindemann, M. J. Manfra, C. Schrade, and C. Schönenberger, Gate-tunable Josephson diode in proximitized InAs supercurrent interferometers, Phys. Rev. Res. 5, 033131 (2023).
  * Roy and Saleur [2019] A. Roy and H. Saleur, Quantum electronic circuit simulation of generalized sine-Gordon models, Phys. Rev. B 100, 155425 (2019).
  * Roy  _et al._ [2021] A. Roy, D. Schuricht, J. Hauschild, F. Pollmann, and H. Saleur, The quantum sine-Gordon model with quantum circuits, Nucl. Phys. B 968, 115445 (2021).
  * Roy and Lukyanov [2023] A. Roy and S. L. Lukyanov, Soliton Confinement in a Quantum Circuit, arXiv (2023), 2302.06289 .
  * Di Francesco  _et al._ [1997] P. Di Francesco, P. Mathieu, and D. Senechal,  _Conformal Field Theory_, Graduate Texts in Contemporary Physics (Springer-Verlag, New York, 1997).
  * Roy [2023] A. Roy, Quantum Electronic Circuits for Multicritical Ising Models, arXiv (2023), 2306.04346 .
  * Mong  _et al._ [2014a] R. S. K. Mong, D. J. Clarke, J. Alicea, N. H. Lindner, P. Fendley, C. Nayak, Y. Oreg, A. Stern, E. Berg, K. Shtengel, and M. P. A. Fisher, Universal Topological Quantum Computation from a Superconductor-Abelian Quantum Hall Heterostructure, Phys. Rev. X 4, 011036 (2014a).
  * Stoudenmire  _et al._ [2015] E. M. Stoudenmire, D. J. Clarke, R. S. K. Mong, and J. Alicea, Assembling Fibonacci anyons from a ${\mathbb{Z}}_{3}$ parafermion lattice model, Phys. Rev. B 91, 235112 (2015).
  * Hu and Kane [2018] Y. Hu and C. L. Kane, Fibonacci topological superconductor, Phys. Rev. Lett. 120, 066801 (2018).
  * Li  _et al._ [2020] C. Li, H. Ebisu, S. Sahoo, Y. Oreg, and M. Franz, Coupled wire construction of a topological phase with chiral tricritical Ising edge modes, Phys. Rev. B 102, 165123 (2020).
  * Friedan  _et al._ [1984] D. Friedan, Z. Qiu, and S. Shenker, Conformal Invariance, Unitarity, and Critical Exponents in Two Dimensions, Phys. Rev. Lett. 52, 1575 (1984).
  * Rahmani  _et al._ [2015a] A. Rahmani, X. Zhu, M. Franz, and I. Affleck, Emergent Supersymmetry from Strongly Interacting Majorana Zero Modes, Phys. Rev. Lett. 115, 166401 (2015a).
  * O’Brien and Fendley [2018] E. O’Brien and P. Fendley, Lattice Supersymmetry and Order-Disorder Coexistence in the Tricritical Ising Model, Phys. Rev. Lett. 120, 206403 (2018).
  * Rahmani  _et al._ [2015b] A. Rahmani, X. Zhu, M. Franz, and I. Affleck, Phase diagram of the interacting Majorana chain model, Phys. Rev. B 92, 235123 (2015b).
  * Zhu and Franz [2016] X. Zhu and M. Franz, Tricritical Ising phase transition in a two-ladder Majorana fermion lattice, Phys. Rev. B 93, 195118 (2016).
  * Ebisu  _et al._ [2019] H. Ebisu, E. Sagi, and Y. Oreg, Supersymmetry in the Insulating Phase of a Chain of Majorana Cooper Pair Boxes, Phys. Rev. Lett. 123, 026401 (2019).
  * Ejima  _et al._ [2016] S. Ejima, F. H. L. Essler, F. Lange, and H. Fehske, Ising tricriticality in the extended Hubbard model with bond dimerization, Phys. Rev. B 93, 235118 (2016).
  * Ejima  _et al._ [2018a] S. Ejima, F. Lange, F. H. Essler, and H. Fehske, Critical behavior of the extended Hubbard model with bond dimerization, Phys. B: Cond. Mat. 536, 474 (2018a).
  * Mong  _et al._ [2014b] R. S. K. Mong, D. J. Clarke, J. Alicea, N. H. Lindner, and P. Fendley, Parafermionic conformal field theory on the lattice, J. Phys. A: Math. Theor. 47, 452001 (2014b).
  * Slagle  _et al._ [2021] K. Slagle, D. Aasen, H. Pichler, R. S. K. Mong, P. Fendley, X. Chen, M. Endres, and J. Alicea, Microscopic characterization of Ising conformal field theory in Rydberg chains, Phys. Rev. B 104, 235109 (2021).
  * Slagle  _et al._ [2022] K. Slagle, Y. Liu, D. Aasen, H. Pichler, R. S. K. Mong, X. Chen, M. Endres, and J. Alicea, Quantum spin liquids bootstrapped from ising criticality in rydberg arrays, Phys. Rev. B 106, 115122 (2022).
  * Buyskikh  _et al._ [2019] A. S. Buyskikh, L. Tagliacozzo, D. Schuricht, C. A. Hooley, D. Pekker, and A. J. Daley, Spin Models, Dynamics, and Criticality with Atoms in Tilted Optical Superlattices, Phys. Rev. Lett. 123, 090401 (2019).
  * Giamarchi [2003] T. Giamarchi,  _Quantum Physics in One Dimension_ (Clarendon Press, Oxford, England, UK, 2003).
  * Haegeman  _et al._ [2011] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Time-dependent variational principle for quantum lattices, Phys. Rev. Lett. 107, 070601 (2011).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94, 165116 (2016).
  * Zauner-Stauber  _et al._ [2018] V. Zauner-Stauber, L. Vanderstraeten, M. T. Fishman, F. Verstraete, and J. Haegeman, Variational optimization algorithms for uniform matrix product states, Phys. Rev. B 97, 045145 (2018).
  * Beenakker [1991] C. W. J. Beenakker, Universal limit of critical-current fluctuations in mesoscopic Josephson junctions, Phys. Rev. Lett. 67, 3836 (1991).
  * Kringhøj  _et al._ [2018] A. Kringhøj, L. Casparis, M. Hell, T. W. Larsen, F. Kuemmeth, M. Leijnse, K. Flensberg, P. Krogstrup, J. Nygård, K. D. Petersson, and C. M. Marcus, Anharmonicity of a superconducting qubit with a few-mode Josephson junction, Phys. Rev. B 97, 060508 (2018).
  * [39] See supplemental material for further data and technical details.
  * Heikkilä  _et al._ [2002] T. T. Heikkilä, J. Särkkä, and F. K. Wilhelm, Supercurrent-carrying density of states in diffusive mesoscopic josephson weak links, Phys. Rev. B 66, 184513 (2002).
  * Zamolodchikov [1986] A. B. Zamolodchikov, Conformal Symmetry and Multicritical Points in Two-Dimensional Quantum Field Theory. (In Russian), Sov. J. Nucl. Phys. 44, 529 (1986).
  * Lässig  _et al._ [1991] M. Lässig, G. Mussardo, and J. L. Cardy, The scaling region of the tricritical Ising model in two dimensions, Nucl. Phys. B 348, 591 (1991).
  * Lepori  _et al._ [2008] L. Lepori, G. Mussardo, and G. Z. Tóth, The particle spectrum of the tricritical Ising model with spin reversal symmetric, J. Stat. Mech.: Theory Exp. 2008 (09), P09004.
  * Tóth [2004] G. Z. Tóth, A nonperturbative study of phase transitions in the multi-frequency sine-Gordon model, J. Phys. A: Math. Gen. 37, 9631 (2004).
  * Orignac and Giamarchi [2001] E. Orignac and T. Giamarchi, Meissner effect in a bosonic ladder, Phys. Rev. B 64, 144515 (2001).
  * Petrescu and Le Hur [2013] A. Petrescu and K. Le Hur, Bosonic Mott Insulator with Meissner Currents, Phys. Rev. Lett. 111, 150601 (2013).
  * Piraud  _et al._ [2015] M. Piraud, F. Heidrich-Meisner, I. P. McCulloch, S. Greschner, T. Vekua, and U. Schollwöck, Vortex and Meissner phases of strongly interacting bosons on a two-leg ladder, Phys. Rev. B 91, 140406 (2015).
  * Petrescu and Le Hur [2015] A. Petrescu and K. Le Hur, Chiral Mott insulators, Meissner effect, and Laughlin states in quantum ladders, Phys. Rev. B 91, 054520 (2015).
  * Greschner  _et al._ [2015] S. Greschner, M. Piraud, F. Heidrich-Meisner, I. P. McCulloch, U. Schollwöck, and T. Vekua, Spontaneous increase of magnetic flux and chiral-current reversal in bosonic ladders: Swimming against the tide, Phys. Rev. Lett. 115, 190402 (2015).
  * Greschner  _et al._ [2016] S. Greschner, M. Piraud, F. Heidrich-Meisner, I. P. McCulloch, U. Schollwöck, and T. Vekua, Symmetry-broken states in a system of interacting bosons on a two-leg ladder with a uniform abelian gauge field, Phys. Rev. A 94, 063628 (2016).
  * Haller  _et al._ [2020] A. Haller, A. S. Matsoukas-Roubeas, Y. Pan, M. Rizzi, and M. Burrello, Exploring helical phases of matter in bosonic ladders, Phys. Rev. Res. 2, 043433 (2020).
  * Fazio and van der Zant [2001] R. Fazio and H. van der Zant, Quantum phase transitions and vortex dynamics in superconducting networks, Phys. Rep. 355, 235 (2001).
  * Goldstein  _et al._ [2013] M. Goldstein, M. H. Devoret, M. Houzet, and L. I. Glazman, Inelastic Microwave Photon Scattering off a Quantum Impurity in a Josephson-Junction Array, Phys. Rev. Lett. 110, 017002 (2013).
  * Giamarchi and Schulz [1987] T. Giamarchi and H. J. Schulz, Localization and Interaction in One-Dimensional Quantum Fluids, Europhys. Lett. 3, 1287 (1987).
  * Orignac and Giamarchi [1998] E. Orignac and T. Giamarchi, Weakly disordered spin ladders, Phys. Rev. B 57, 5812 (1998).
  * Crépin  _et al._ [2011] F. Crépin, N. Laflorencie, G. Roux, and P. Simon, Phase diagram of hard-core bosons on clean and disordered two-leg ladders: Mott insulator–Luttinger liquid–Bose glass, Phys. Rev. B 84, 054517 (2011).
  * Delfino and Mussardo [1998] G. Delfino and G. Mussardo, Non-integrable aspects of the multi-frequency sine-Gordon model, Nucl. Phys. B 516, 675 (1998).
  * Tagliacozzo  _et al._ [2008] L. Tagliacozzo, Thiago. R. de Oliveira, S. Iblisdir, and J. I. Latorre, Scaling of entanglement support for matrix product states, Phys. Rev. B 78, 024410 (2008).
  * Zanardi and Paunković [2006] P. Zanardi and N. Paunković, Ground state overlap and quantum phase transitions, Phys. Rev. E 74, 031123 (2006).
  * Cozzini  _et al._ [2007] M. Cozzini, R. Ionicioiu, and P. Zanardi, Quantum fidelity and quantum phase transitions in matrix product states, Phys. Rev. B 76, 104420 (2007).
  * Rams and Damski [2011] M. M. Rams and B. Damski, Quantum fidelity in the thermodynamic limit, Phys. Rev. Lett. 106, 055701 (2011).
  * Rossini and Vicari [2018] D. Rossini and E. Vicari, Ground-state fidelity at first-order quantum transitions, Phys. Rev. E 98, 062137 (2018).
  * Kane and Fisher [1996] C. L. Kane and M. P. A. Fisher, Thermal Transport in a Luttinger Liquid, Phys. Rev. Lett. 76, 3192 (1996).
  * Bernard and Doyon [2016] D. Bernard and B. Doyon, Conformal field theory out of equilibrium: a review, J. Stat. Mech.: Theory Exp. 2016 (6), 064005.
  * Gawędzki  _et al._ [2018] K. Gawędzki, E. Langmann, and P. Moosavi, Finite-Time Universality in Nonequilibrium CFT, J. Stat. Phys. 172, 353 (2018).
  * Ejima  _et al._ [2018b] S. Ejima, T. Yamaguchi, F. H. L. Essler, F. Lange, Y. Ohta, and H. Fehske, Exotic criticality in the dimerized spin-1 $XXZ$ chain with single-ion anisotropy, SciPost Phys. 5, 059 (2018b).
  * Feiguin  _et al._ [2007] A. Feiguin, S. Trebst, A. W. W. Ludwig, M. Troyer, A. Kitaev, Z. Wang, and M. H. Freedman, Interacting Anyons in Topological Quantum Liquids: The Golden Chain, Phys. Rev. Lett. 98, 160409 (2007).
  * Zou  _et al._ [2018] Y. Zou, A. Milsted, and G. Vidal, Conformal data and renormalization group flow in critical quantum spin chains using periodic uniform matrix product states, Phys. Rev. Lett. 121, 230402 (2018).
  * Cubero  _et al._ [2022] A. C. Cubero, R. M. Konik, M. Lencsés, G. Mussardo, and G. Takács, Duality and form factors in the thermally deformed two-dimensional tricritical Ising model, SciPost Phys. 12, 162 (2022).
  * Lencsés  _et al._ [2022] M. Lencsés, G. Mussardo, and G. Takács, Confinement in the tricritical Ising model, Phys. Lett. B 828, 137008 (2022).
  * Centre [2021] J. S. Centre, Jureca: Data centric and booster modules implementing the modular supercomputing architecture at jülich supercomputing centre, Journal of large-scale research facilities 7, A182 (2021).
  * Beenakker [1997] C. W. J. Beenakker, Random-matrix theory of quantum transport, Rev. Mod. Phys. 69, 731 (1997).
  * Nichele  _et al._ [2020] F. Nichele, E. Portolés, A. Fornieri, A. M. Whiticar, A. C. C. Drachmann, S. Gronin, T. Wang, G. C. Gardner, C. Thomas, A. T. Hatke, M. J. Manfra, and C. M. Marcus, Relating Andreev Bound States and Supercurrents in Hybrid Josephson Junctions, Phys. Rev. Lett. 124, 226801 (2020).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Mussardo [2020] G. Mussardo,  _Statistical Field Theory_ (Oxford University Press, Oxford, England, UK, 2020).
  * Calabrese and Cardy [2004] P. Calabrese and J. Cardy, Entanglement entropy and quantum field theory, J. Stat. Mech.: Theory Exp. 2004 (06), P06002.
  * Fishman  _et al._ [2022] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022).