††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Superconducting qubit based on twisted cuprate van der Waals
heterostructures

Valentina Brosco  Institute for Complex Systems (ISC) Consiglio Nazionale
delle Ricerche and Physics Department University of Rome, “La Sapienza”, P.le
A. Moro, 2 (00185) Roma, Italy  Giuseppe Serpico  Max Planck Institute for
Chemical Physics of Solids, 01187 Dresden, Germany  Department of Physics,
University of Naples Federico II, Via Cintia, Naples, 80126, Italy  Valerii
Vinokur  Terra Quantum AG, Kornhausstrasse 25, CH-9000 St. Gallen, Switzerland
Physics Department, CUNY, City College of City University of New York, 160
Convent Ave, New York, NY 10031, USA  Nicola Poccia  Leibniz Institute for
Solid State and Materials Science Dresden (IFW Dresden), 01069 Dresden,
Germany  Uri Vool  Max Planck Institute for Chemical Physics of Solids, 01187
Dresden, Germany

###### Abstract

Van-der-Waals (vdW) assembly enables the fabrication of novel Josephson
junctions featuring an atomically sharp interface between two exfoliated and
relatively twisted $\rm{Bi_{2}Sr_{2}CaCu_{2}O_{8+x}}$ (Bi2212) flakes. In a
range of twist angles around $45^{\circ}$, the junction provides a regime
where the interlayer two-Cooper pair tunneling dominates the current-phase
relation. Here we propose employing this novel junction to realize a
capacitively shunted qubit that we call flowermon. The $d$-wave nature of the
order parameter endows the flowermon with inherent protection against charge-
noise-induced relaxation and quasiparticle-induced dissipation. This
inherently protected qubit paves the way to a new class of high-coherence
hybrid superconducting quantum devices based on unconventional
superconductors.

Superconducting microwave circuits are macroscopic devices that mimic the
quantum properties of atoms. Thanks to their flexible design and strong
coupling, superconducting circuits lead the implementation of hardware for
quantum technology Devoret and Schoelkopf (2013). Currently, the most commonly
used superconducting circuits are based on the transmon qubit, having a simple
and robust design with a single Josephson junction and a shunting capacitor
Koch et al. (2007). However, the plasmonic nature of the transmon and the
macroscopic size of the shunt capacitor plates reduce its coherence time and
enhance cross-talk between qubits, posing significant limitations for the
implementation of advanced quantum devices comprising a large number of qubits
and requiring low error rates Brink et al. (2018); Arute et al. (2019). As a
result, a growing number of works explores alternative designs with inherent
protection, such as the rhombus Blatter et al. (2001); Douçot et al. (2003);
Gladchenko et al. (2009); Bell et al. (2014), fluxonium Manucharyan et al.
(2009); Nguyen et al. (2019), bifluxonKalashnikov et al. (2020), blochnium
Pechenezhskiy et al. (2020), KITE Smith et al. (2020, 2022), $0-\pi$ Brooks et
al. (2013); Groszkowski et al. (2018); Gyenis et al. (2021) and nanowire-based
Larsen et al. (2020); Ciaccia et al. (2023) qubits. While significantly higher
coherence times are predicted for such circuits in the ideal case, their
implementation relies on multiple junctions in a flux loop, making them
vulnerable to flux noise and inevitable fabrication imperfections. Here we
develop a novel single-junction qubit with inherent protection that avoids
these drawbacks and provides a robust platform for future quantum devices.

Figure 1: Potential design for the flowermon qubit. (a) A relative twist of
two d-wave flakes placed together to form a Josephson junction can suppress
Cooper pair tunneling due to momentum mismatch. At $45^{\circ}$ the mismatch
completely suppresses single Cooper pair tunneling, and two-pair tunneling
dominates the junction. (b) The design of the flowermon with a single $d$-wave
junction shunted by a large capacitor similar to the transmon qubit. The 3D
design shows a possible physical implementation with the capacitor pads of a
conventional superconductor coupled to the junction. (c) Josephson’s potential
for different values of the twist angle.

Pioneering works Ioffe et al. (1999); Blais and Zagoskin (2000); Blatter et
al. (2001); Bauch et al. (2006) proposed to utilize the suppression of
tunneling in the $d$-wave based Josephson junctions to realize a qubit with an
enhanced coherence. Note that since a single Cooper-pair tunneling between
twisted superconducting $d$-wave islands having different orientations is
suppressed due to the momentum mismatch, see Fig. 1a, the contribution of the
second-harmonic $\cos(2\varphi)$ of the Josephson energy becomes dominant, and
the consequent degeneracy of the spectrum can be exploited for creating a
decoherence-protected qubit. In the standard approaches using the grain
orientation mismatch, the Josephson junctions at grain boundaries of
YBa2Cu3O7-x with relative granule mismatch angle close to $\pi/4$ demonstrated
the degeneracy of the ground-state Ilichev et al. (2001). However, fabrication
complexity and the resulting low junction quality hindered the implementation
of quantum devices based on grain boundary junctions.

Recent technological progress enabled to preserve nearly perfect
superconductivity Yu et al. (2019); Zhao et al. (2019) and lattice structure
Poccia et al. (2020) in the isolated atomically thin Bi2Sr2CaCu2O8+x (Bi2212)
crystals. This enabled creation of the Bi2212-based junctions, employing a
novel technique that allowed to control the diffusion of oxygen interstitials,
which is the main source of detrimental disorder in cuprate Bi2212. A
cryogenic stacking protocol has been developed utilizing the freezing of the
motion of the oxygen interstitials at temperatures below 200 K Poccia et al.
(2011), which allowed manufacturing the vdW flakes and creation of an
atomically sharp interface between Bi2212 crystals. Such protocol enables the
fabrication of high-quality Josephson junctions, displaying a strong
dependence of the Josephson energy on the twist angle Zhao et al. (2021); Lee
et al. (2023); Martini et al. (2023).

Here we propose to utilize such a junction as a platform for a novel highly-
protected qubit that we call flowermon. We focus on the simplest design
comprising a single capacitively shunted twisted vdW Josephson junction, as
sketched in Fig. 1(b). We quantitatively analyze the behavior of the flowermon
qubit, highlighting possible decoherence mechanisms for experimentally-
available parameters, and discuss possible control and measurement schemes.

The flowermon junction has a twist angle $\theta$ between the axis associated
with the $d$-wave order parameters of two overlapping superconducting flakes,
see Fig. 1a, and is shunted by a capacitor with charging energy
$E_{\mathrm{C}}=e^{2}/(2C)$ where $C$ is the shunting capacitance which
dominates over the internal junction capacitance, see Fig. 1b.

Figure 2: Flowermon low-energy spectrum. (a) Energy levels of the flowermon
qubit as a function of the angle $\theta$. Even and odd levels are shown in
red and blue lines respectively; as $\theta$ approaches $\pi/4$, the levels
combine into quasi-degenerate doublets. Here we assume
$E_{\mathrm{\kappa}}/E_{\text{J}}=0.1$ and $E_{J}/E_{\text{C}}=2000$. (b) and
(c) Level structure and potential energy for $\theta=30^{\circ}$ and
$\theta=43^{\circ}$. (d) and (e) The wave-functions for the ground and first
excited state $|\Psi_{0}\rangle$ and $|\Psi_{1}\rangle$ in the phase basis for
$\theta=30^{\circ}$ and $\theta=43^{\circ}$. (f) and (g) Evolution of the
structure of $|\Psi_{0}\rangle$ and $|\Psi_{1}\rangle$ in the charge basis as
a function of the twisting angle $\theta$.

The circuit Hamiltonian is

$\displaystyle
H=4E_{\text{C}}(\hat{n}-n_{\text{g}})^{2}-E_{\text{J}\theta}\\!\cos(\hat{\varphi})+E_{\mathrm{\kappa}}\\!\cos(2\hat{\varphi})\,,$
(1)

where $\hat{\varphi}$ denotes the phase difference across the junction and
$\hat{n}$ is the conjugate charge. In the above equation
$E_{\mathrm{J}\theta}$ and $E_{\mathrm{\kappa}}$ represent, respectively, the
energy associated with the coherent tunneling of single and double Cooper
pairs across the junction, while the gate charge $n_{\text{g}}$ accounts for
charge fluctuations induced by external electric fields. The competition
between the two Josephson tunneling terms arises due to the peculiar
dependence of the Josephson energy on the twisting angle,
$E_{\mathrm{J}\theta}=E_{\mathrm{J}}\cos(2\theta)$, associated with the
$d$-wave structure of the order parameter as has been recently highlighted in
Refs. Can et al. (2021); Tummuru et al. (2022) and confirmed experimentally in
Ref. Zhao et al. (2021). The shape of the Josephson potential,
$U_{\mathrm{J}}(\varphi)$ as a function of the twisting angle, is shown in
Fig. 1(c). As $\theta$ increases from $0$ to $\pi/4$, the Josephson potential
develops a symmetric double-well structure with two minima at
$\varphi=\pm\varphi_{0}$ with
$\varphi_{0}=\arccos(E_{\mathrm{J}\theta}/4E_{\mathrm{\kappa}})$ separated by
a barrier

$\Delta
U_{\mathrm{J}}=\left(E_{\mathrm{J}\theta}-4E_{\kappa}\right)^{2}/8E_{\kappa}.$
(2)

Using Eq. (2), we express the critical angle at which the double-well
structure arises as
$\theta_{\mathrm{c}}=\frac{1}{2}\arccos(4E_{\kappa}/E_{\mathrm{J}})$. The
ratio $E_{\mathrm{\kappa}}/E_{\mathrm{J}}$ thus determines the range of angles
in which the qubit can be realized. There is some discrepancy between the
theoretical estimates of this ratio Can et al. (2021) as 0.1-0.2, and
experimental observations Zhao et al. (2021) giving 0.02-0.05. Throughout the
paper, we use $E_{\kappa}/E_{\mathrm{J}}=0.1$; other values will lead to the
renormalization of the operating angles. Note that
$E_{\mathrm{J}}\propto\Delta$ and $E_{\mathrm{\kappa}}\propto\Delta^{2}$,
hence high-temperature superconductors such as Bi2212 have a natural advantage
due to higher two-Copper-pair tunneling rates. Additional effects such as
inhomogeneity in the twisted junction can further increase
$E_{\mathrm{\kappa}}$ Yuan et al. (2023). For further discussion of the
relation between junction and qubit parameters, see supplementary information
sup .

For $\theta_{\mathrm{c}}<\theta\lesssim\pi/4$ and sufficiently small
$E_{\mathrm{C}}$, the low-energy spectrum reduces to a set of quasi-degenerate
doublets, see Fig. 2(a). Splitting of each doublet is associated with the
tunneling of the phase between the two minima and thus scales exponentially
with the ratio $\Delta U_{\mathrm{J}}/E_{\mathrm{C}}$ Smith et al. (2020); sup
. This transition is shown in Fig. 2(b),(c) where at $\theta=30^{\circ}$ the
spectrum is quasi-harmonic, reminiscent of the spectrum of the transmon qubit,
while at $\theta=43^{\circ}$ the double-well fully develops, and the lowest
energy levels are doubly degenerate. In the latter case, the phase is not able
to tunnel between two adjacent minima, and, as shown in Fig. 2(e), the
corresponding wave functions are localized in two energy wells. The quasi-
degeneracy of the levels is associated with the symmetry of the wave
functions. As the angle approaches $\pi/4$ the wave functions of the even and
odd energy levels develop a distinct parity in the charge basis. Specifically,
at high twisting angles, $\theta\gtrsim 40^{\circ}$, the ground state
$|\psi_{0}\rangle$ contains only even Cooper pair number states, while the
first excited state $|\psi_{1}\rangle$ contains only odd Cooper pair number
states, see Fig. 2(f),(g). Therefore, when the Josephson energy is completely
dominated by double Cooper pair tunneling, the Cooper pair number parity is
conserved.

To illustrate the relevance of the flowermon symmetry properties for qubit
protection, we consider the expression of the relaxation rate induced by the
capacitive losses Makhlin et al. (2001)

$\Gamma_{1c}=\frac{(8E_{\mathrm{C}})^{2}}{\hbar^{2}}S_{n_{g}}(\omega_{01})n_{01}^{2}\,,$
(3)

where $n_{xy}=\langle\psi_{x}|\hat{n}|\psi_{y}\rangle$ is the matrix element
of the charge operator and $S_{n_{g}}(\omega_{01})$ is the spectral density of
the capacitive noise at the qubit operational frequency. Due to the very small
overlap in the charge basis $\hat{n}$ between $|\psi_{0}\rangle$ and
$|\psi_{1}\rangle$ in the flowermon limit, the corresponding matrix element is
exponentially suppressed, preventing energy relaxation, see Fig. 3(b). Looking
at twisting angles $\theta\gtrsim 40^{\circ}$, we, thus, expect an enhancement
of the relaxation time over several orders of magnitude as compared to the
standard transmon devices which are typically limited by such dielectric
losses Wang et al. (2015); Ahmad et al. (2022).

The corresponding dephasing rate due to charge fluctuations can be estimated
as

$\Gamma_{\mathrm{\varphi}c}\simeq\frac{32E_{C}^{2}}{\hbar^{2}}S_{n_{g}}(0)|n_{11}-n_{00}|^{2}\,,$
(4)

where $S_{n_{g}}(0)$ is the spectral density of the charge noise at zero
frequency. In the case where $E_{\mathrm{C}}$ is sufficiently small with
respect to $E_{\mathrm{J}}$ and $E_{\mathrm{\kappa}}$, the coefficient
$n_{11}-n_{00}$ is exponentially suppressed. The suppression factor can be
quantitatively estimated by mapping the Hamiltonian onto a tunneling problem,
see Supplementary Information sup . Thanks to the large shunt capacitor, the
flowermon inherits protection against dephasing induced by the charge noise
from the transmon.

Figure 3: Protection from decoherence. (a) The charge matrix elements
$n_{01}^{2}$ and $|n_{11}-n_{00}|^{2}$, related to charge/capacitive
relaxation and dephasing rates respectively, plotted vs.
$E_{\mathrm{J}}/E_{\mathrm{C}}$ and the twisting angle $\theta$. Different
regimes are shown, with the flowermon exhibiting protection from both
relaxation and dephasing processes. (b) Charge matrix elements $n_{01}$,
$n_{03}$, $n_{23}$, and $n_{12}$ vs. the twisting angle $\theta$ and
$E_{\mathrm{J}}/E_{\mathrm{C}}=2000$. Both $n_{01}$ and $n_{23}$ decay
exponentially as $\theta$$\rightarrow$$45^{\circ}$ but with a significant
offset. (c) Dependence of the quasiparticle relaxation rate on the ratio
$\Delta/k_{\mathrm{B}}T$ for different values of the twisting angles ranging
between $\theta=0^{\circ}$ and $\theta=40^{\circ}$.

.

In superconducting qubits based on the $s$-wave junctions, the existence of an
energy gap guarantees that the decoherence induced by thermally-equilibrium
quasiparticles is exponentially suppressed Catelani et al. (2011). Here we
show that for sufficiently clean interfaces, an analogous gap behavior is also
present in the flowermon at high twisting angles. To that end, we write the
junction Hamiltonian as Can et al. (2021)

$H=H_{\mathrm{L}}+H_{\mathrm{R}}+H_{\mathrm{T}}\,,$ (5)

where $H_{\mathrm{L}}$ and $H_{\mathrm{R}}$ represent the BCS Hamiltonians of
the two layers, i.e.,

$\displaystyle H_{\mathrm{L,R}}=\sum_{{\bf k}\sigma}\xi_{{\bf
k}{\mathrm{L,R}}}c^{\dagger}_{{\bf k}\sigma\mathrm{L,R}}c_{{\bf
k}\sigma{\mathrm{L,R}}}+$ $\displaystyle+\sum_{\bf k}\left(\Delta_{{\bf
k}\mathrm{L,R}}c^{\dagger}_{{\bf
k}\uparrow\mathrm{L,R}}c^{\dagger}_{{\bf-k}\downarrow\mathrm{L,R}}+{\rm
H.c.}\right)\,,$ (6)

with $c^{\dagger}_{{\bf k}\sigma\mathrm{L,R}}$ and $c_{{\bf
k}\sigma\mathrm{L,R}}$ denoting creation and annihilation operators of
electrons with spin $\sigma$ and momentum ${\bf k}$ in layers $L$ and $R$,
while $H_{\mathrm{T}}$ represents the interlayer tunneling

$H_{\mathrm{T}}=\sum_{{\bf kp}\sigma}\left(t_{\bf
kp}e^{-i\varphi/2}c^{\dagger}_{{\bf k}\sigma{\mathrm{L}}}c_{{\bf
p}\sigma{\mathrm{R}}}+{\rm H.c.}\right).$ (7)

The Hamiltonian $H$ thus depends on a small number of parameters, namely, the
$d$-wave gap functions defined as

$\Delta_{{\bf k}{\mathrm{L}}}=\Delta_{d}\cos(2\theta_{\bf k})\,\,{\rm
and}\,\,\Delta_{{\bf k}{\mathrm{R}}}=\Delta_{d}\cos(2\theta_{\bf k}-2\theta),$
(8)

with $\theta_{\bf k}$ denoting the polar angle in the plane $k_{x},k_{y}$, the
electronic dispersion around the Fermi surface $\xi_{\bf
k}=\hbar^{2}k^{2}/(2m)-\mu$, and the interlayer tunnel amplitude $t_{\bf kp}$.
The latter plays a central role in determining the junction’s properties, as
discussed in Ref. Haenel et al. (2022). Following Refs. Yokoyama et al.
(2007); Can et al. (2021); Bruder et al. (1995), we consider the coherent
tunneling limit with the in-plane momentum conservation, $t_{\bf
kp}=t\,\delta_{\bf k_{\parallel}p_{\parallel}}$, suitable to describe weakly
disordered $c$-axis junctions Bruder et al. (1995). More accurate microscopic
treatment of $t_{\bf kp}$, including momentum dependence, has been previously
considered Barash et al. (1995) and can affect the Josephson coupling terms
themselves in addition to quasiparticle tunneling Song et al. (2022). This
full treatment is beyond the scope of this work and will be the subject of a
forthcoming publication.

We calculate the quasiparticle-induced relaxation rate using the approach
developed by Refs. Catelani et al. (2011, 2012); Catelani (2014), which can be
expressed as sup

$\Gamma_{\downarrow,{\rm
qp}}=t^{2}|\langle\psi_{0}|\sin(\varphi/2)|\psi_{1}\rangle|^{2}S^{\rm\theta}_{\rm
qp}(\omega_{01}),$ (9)

where $S^{\rm\theta}_{\rm qp}(\omega_{01})$ is the spectral density of the
noise fluctuations induced by quasiparticle tunneling. In the above equation,
the $\theta$ dependence of the relaxation rate can be traced back to the
dependence of both the matrix element
$\langle\psi_{0}|\sin(\varphi/2)|\psi_{1}\rangle$ and spectral density on the
twisting angle.

Figure 4: Manipulation of the flowermon logical states. (a) Logical states
manipulation can be performed via the 2-nd and 3-rd excited states in a regime
where the 0-1 matrix element and frequency are near 0, but the amplitude of
the 2$\leftrightarrow$3 transition is still finite. The control is performed
by standard microwave $\pi$-pulses between the 1-2, and 0-3 states, the
desired control gate $U$ pulse between the 2-3 states, and again $\pi$-pulses
between the 1-2, and 0-3 states. Alternatively, the 0-1 transition can be
driven indirectly by a Raman process via simultaneous off-resonant driving of
the 0-3, 2-3, and 2-1 transitions (orange arrows). Note that this protocol
works even when $\omega_{01}\rightarrow 0$ because the selection rule forbids
0-2 and 1-3 transitions. (b) A sketch of a possible implementation of the
device with control and readout lines.

.

Figure 3(c) shows the results of the numerical calculation of the
quasiparticle relaxation rate given by Eq. (9). At high temperatures, high
twisting angles correspond to a larger quasiparticle relaxation rate due to
the larger prefactor $\langle\psi_{0}|\sin(\varphi/2)|\psi_{1}\rangle$. On the
contrary, in the low-temperature limit, high twisting angles yield a strong
exponential suppression of the quasi-particle rate $\Gamma_{\downarrow,{\rm
qp}}$, as nodal quasiparticles are forbidden from tunneling due to a momentum
mismatch. The numerical results are well fitted by the simple analytical
formula, $S^{\rm\theta}_{\rm
qp}(\omega_{01})\propto\,e^{-\Delta_{d}/(k_{B}T)\sin(2\theta)^{2}}$, where, as
mentioned above, $\Delta_{d}$ is the $d$-wave gap in the quasiparticle
spectrum at zero twisting angle; this result is reminiscent of what was found
in Ref. Yokoyama et al. (2007), see supplementary material sup  for details.
This fit indicates that under appropriate conditions the quasi-particles in
twisted cuprate vdW heterostructures are effectively gapped, in similarity to
what happens in the s-wave junctions.

So far we considered only thermal quasiparticles but superconducting devices
are often limited by non-equilibrium quasiparticles Martinis et al. (2009); de
Visser et al. (2014); Vool et al. (2014); Wang et al. (2014); Serniak et al.
(2018); Vepsäläinen et al. (2020); Cardani et al. (2021). We expect that the
gap behavior highlighted here will similarly protect the flowermon from the
decoherence induced by non-equilibrium quasiparticle tunneling Pop et al.
(2014). Note that similar protection can be obtained by an s/d-wave Josephson
junction Ioffe et al. (1999); Patel et al. (2023).

The decoupling of the flowermon qubit from noise sources hampers the direct
control and readout of the qubit since the qubit energy $\hbar\omega_{01}$
(see Fig. 2(a)) and the charge matrix element $n_{01}$ (Fig. 3(b)) are
exponentially suppressed as $\theta$$\rightarrow$$45^{\circ}$. However, the
structure of the flowermon spectrum allows for manipulation through the higher
energy levels without sacrificing protection, see Fig. 4(a). Specifically,
there is a wide range of angles where the qubit matrix element $n_{01}$ is
suppressed while the coupling between the 2-nd and 3-rd excited states
$n_{23}$ is still finite (Fig. 3b), e.g., $\theta=40^{\circ}$ corresponds to
$n_{01}=0.01$ and $n_{23}=0.4$. Thus, one can use the excited states for
manipulation. For example, a $\pi/2$-pulse in the logical 0-1 space can be
performed using standard microwave $\pi$-pulses between the 1-2, and 0-3
states, a $\pi/2$-pulse between the 2-3 states, and again $\pi$-pulses between
the 1-2, and 0-3 states. To avoid populating the higher states, 0-1
manipulation can also be performed via a simultaneous application of the off-
resonant Raman drives through the same energy levels Vool et al. (2018), see
the orange line in Fig. 4(a). Note that this scheme works even if
$\omega_{01}$ is smaller than the frequency selectivity of the control pulses
since the transitions 0-2, 1-3 are forbidden by the selection rules
$n_{02}=n_{13}=0$. The higher excited states can be used similarly for
measuring the state of the qubit, either by measuring non-protected
transitions such as 0-3, 1-2 directly Cottet et al. (2021) or through the
dispersive coupling of these transitions to a cavity mode Smith et al. (2020).
Figure 4(b) shows a possible experimental implementation of such a qubit with
the control and readout lines. A more quantitative discussion of this design
with realistic qubit parameters is given in sup .

In conclusion, we have developed a qubit design based on a twisted cuprate
heterostructure, which we call the flowermon. The simple design of the
flowermon, a single junction shunted by a large capacitance similar to that of
a transmon, does not require flux tuning or precise control of the fabrication
parameters. However, the $d$-wave structure of the superconducting phase in
the unique junction endows the qubit with physical protection against
decoherence. We have discussed the low energy spectrum structure of the qubit
at various twist angles and quantitatively analyzed its sensitivity to
decoherence mechanisms. Finally, we have presented schemes for the
manipulation and readout of this qubit. At twisting angles close to
$45^{\circ}$, the flowermon shows exponential suppression of the charge-
induced noise, as well as protection from nodal quasiparticle tunneling. The
flowermon is thus expected to provide the orders of magnitude improvement in
the coherence time, promising remarkable progress in future superconducting
quantum hardware. Furthermore, the protection offered by the qubit is directly
related to the microscopic properties of the cuprate superconductor. Thus, the
flowermon promises to become the prototype for a new class of hybrid devices
which combine the benefits of quantum materials and coherent quantum circuits.

Acknowledgements– The work is partially supported by the Deutsche
Forschungsgemeinschaft (DFG 452128813, DFG 512734967, DFG 492704387, and DFG
460444718) and co-funded by the European Union (ERC, cQEDscope, 101075962).
The work of V.M.V. is supported by Terra Quantum Inc. and partially by the US
NSF Grant Awards No. NSF 1809188 and 2105048. The authors are deeply grateful
to Shu Yang Frank Zhao, Philip Kim, Bernard van Heck, Carlo Di Castro, Rosario
Fazio, Andrei Zaikin, Kornelius Nielsch, Francesco Tafuri, Giampiero Pepe,
Domenico Montemurro, Davide Massarotti, and Pavel A. Volkov for illuminating
discussions.

## References

  * Devoret and Schoelkopf (2013) M. H. Devoret and R. J. Schoelkopf, Science 339, 1169 (2013).
  * Koch et al. (2007) J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 76, 042319 (2007).
  * Brink et al. (2018) M. Brink, J. M. Chow, J. Hertzberg, E. Magesan, and S. Rosenblatt, 2018 IEEE International Electron Devices Meeting (IEDM) pp. 6.1.1–6.1.3 (2018).
  * 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, et al., Nature 574, 505 (2019).
  * Blatter et al. (2001) G. Blatter, V. B. Geshkenbein, and L. B. Ioffe, Physical Review B 63, 174511 (2001).
  * Douçot et al. (2003) B. Douçot, M. V. Feigel’man, and L. B. Ioffe, Physical Review Letters 90 (2003).
  * Gladchenko et al. (2009) S. Gladchenko, D. Olaya, E. Dupont-Ferrier, B. Douçot, L. B. Ioffe, and M. E. Gershenson, Nature Physics 5 (2009).
  * Bell et al. (2014) M. T. Bell, J. Paramanandam, L. B. Ioffe, and M. E. Gershenson, Physical Review Letters 112, 167001 (2014).
  * Manucharyan et al. (2009) V. E. Manucharyan, J. Koch, L. I. Glazman, and M. H. Devoret, Science 326, 113 (2009).
  * Nguyen et al. (2019) L. B. Nguyen, Y.-H. Lin, A. Somoroff, R. Mencia, N. Grabon, and V. E. Manucharyan, Physical Review X 9, 041041 (2019).
  * Kalashnikov et al. (2020) K. Kalashnikov, W. T. Hsieh, W. Zhang, W.-S. Lu, P. Kamenov, A. Di Paolo, A. Blais, M. E. Gershenson, and M. Bell, PRX Quantum 1, 010307 (2020).
  * Pechenezhskiy et al. (2020) I. V. Pechenezhskiy, R. A. Mencia, L. B. Nguyen, Y.-H. Lin, and V. E. Manucharyan, Nature 585, 368 (2020).
  * Smith et al. (2020) W. C. Smith, A. Kou, X. Xiao, U. Vool, and M. H. Devoret, npj Quantum Information 6 (2020).
  * Smith et al. (2022) W. C. Smith, M. Villiers, A. Marquet, J. Palomo, M. R. Delbecq, T. Kontos, P. Campagne-Ibarcq, B. Douçot, and Z. Leghtas, Physical Review X 12, 021002 (2022).
  * Brooks et al. (2013) P. Brooks, A. Kitaev, and J. Preskill, Physical Review A 87, 052306 (2013).
  * Groszkowski et al. (2018) P. Groszkowski, A. D. Paolo, A. L. Grimsmo, A. Blais, D. I. Schuster, A. A. Houck, and J. Koch, New Journal of Physics 20, 043053 (2018).
  * Gyenis et al. (2021) A. Gyenis, P. S. Mundada, A. Di Paolo, T. M. Hazard, X. You, D. I. Schuster, J. Koch, A. Blais, and A. A. Houck, PRX Quantum 2, 010339 (2021).
  * Larsen et al. (2020) T. W. Larsen, M. E. Gershenson, L. Casparis, A. Kringhøj, N. J. Pearson, R. P. G. McNeil, F. Kuemmeth, P. Krogstrup, K. D. Petersson, and C. M. Marcus, Physical Review Letters 125, 056801 (2020).
  * Ciaccia et al. (2023) C. Ciaccia, R. Haller, A. C. C. Drachmann, T. Lindemann, M. J. Manfra, C. Schrade, and C. Schönenberger, arXiv2306.05467 (2023).
  * Ioffe et al. (1999) L. B. Ioffe, V. B. Geshkenbein, M. V. Feigel’man, A. L. Fauchère, and G. Blatter, Nature 398, 679 (1999).
  * Blais and Zagoskin (2000) A. Blais and A. M. Zagoskin, Physical Review A 61, 042308 (2000).
  * Bauch et al. (2006) T. Bauch, T. Lindström, F. Tafuri, G. Rotoli, P. Delsing, T. Claeson, and F. Lombardi, Science 311, 57 (2006).
  * Ilichev et al. (2001) E. Ilichev, M. Grajcar, R. Hlubina, R. P. J. IJsselsteijn, H. E. Hoenig, H.-G. Meyer, A. Golubov, M. H. S. Amin, A. M. Zagoskin, A. N. Omelyanchouk, et al., Physical Review Letters 86, 5369 (2001).
  * Yu et al. (2019) Y. Yu, L. Ma, P. Cai, R. Zhong, C. Ye, J. Shen, G. D. Gu, X. H. Chen, and Y. Zhang, Nature 575, 156 (2019).
  * Zhao et al. (2019) S. F. Zhao, N. Poccia, M. G. Panetta, C. Yu, J. W. Johnson, H. Yoo, R. Zhong, G. Gu, K. Watanabe, T. Taniguchi, et al., Physical review letters 122, 247001 (2019).
  * Poccia et al. (2020) N. Poccia, S. Y. F. Zhao, H. Yoo, X. Huang, H. Yan, Y. S. Chu, R. Zhong, G. Gu, C. Mazzoli, K. Watanabe, et al., Physical Review Materials 4, 114007 (2020).
  * Poccia et al. (2011) N. Poccia, M. Fratini, A. Ricci, G. Campi, L. Barba, A. Vittorini-Orgeas, G. Bianconi, G. Aeppli, and A. Bianconi, Nature materials 10, 733 (2011).
  * Zhao et al. (2021) S. Y. F. Zhao, N. Poccia, X. Cui, P. A. Volkov, H. Yoo, R. Engelke, Y. Ronen, R. Zhong, G. Gu, S. Plugge, et al., arXiv2108.13455 (2021). URL https://doi.org/10.48550/arXiv.2108.13455.
  * Lee et al. (2023) Y. Lee, M. Martini, T. Confalone, S. Shokri, C. N. Saggau, D. Wolf, G. Gu, K. Watanabe, T. Taniguchi, D. Montemurro, et al., Advanced Materials 35, 2209135 (2023).
  * Martini et al. (2023) M. Martini, Y. Lee, T. Confalone, S. Shokri, C. N. Saggau, D. Wolf, G. Gu, K. Watanabe, T. Taniguchi, D. Montemurro, et al., Materials Today(2023). URL 10.1016/j.mattod.2023.06.007
  * Can et al. (2021) O. Can, T. Tummuru, R. P. Day, I. Elfimov, A. Damascelli, and M. Franz, Nature Physics 17, 519 (2021).
  * Tummuru et al. (2022) T. Tummuru, S. Plugge, and M. Franz, Physical Review B 105, 064501 (2022).
  * Yuan et al. (2023) A. C. Yuan, Y. Vituri, E. Berg, B. Spivak, and S. A. Kivelson, arXiv2305.15472 (2023).
  * (34) See supplementary material.
  * Makhlin et al. (2001) Y. Makhlin, G. Schön, and A. Shnirman, Rev. Mod. Phys. 73, 357 (2001).
  * Wang et al. (2015) C. Wang, C. Axline, Y. Y. Gao, T. Brecht, Y. Chu, L. Frunzio, M. H. Devoret, and R. J. Schoelkopf, Applied Physics Letters 107, 162601 (2015), ISSN 0003-6951.
  * Ahmad et al. (2022) H. G. Ahmad, V. Brosco, A. Miano, L. Di Palma, M. Arzeo, D. Montemurro, P. Lucignano, G. P. Pepe, F. Tafuri, R. Fazio, et al., Phys. Rev. B 105, 214522 (2022).
  * Catelani et al. (2011) G. Catelani, R. J. Schoelkopf, M. H. Devoret, and L. I. Glazman, Physical Review B 84, 064517 (2011).
  * Haenel et al. (2022) R. Haenel, T. Tummuru, and M. Franz, Physical Review B 106, 104505 (2022).
  * Yokoyama et al. (2007) T. Yokoyama, S. Kawabata, T. Kato, and Y. Tanaka, Physical Review B 76, 134501 (2007).
  * Bruder et al. (1995) C. Bruder, A. van Otterlo, and G. T. Zimanyi, Phys. Rev. B 51, 12904 (1995), URL https://link.aps.org/doi/10.1103/PhysRevB.51.12904.
  * Barash et al. (1995) Y. S. Barash, A. V. Galaktionov, and A. D. Zaikin, Phys. Rev. B 52, 665 (1995).
  * Song et al. (2022) X.-Y. Song, Y.-H. Zhang, and A. Vishwanath, Physical Review B 105, L201102 (2022).
  * Catelani et al. (2012) G. Catelani, S. E. Nigg, S. M. Girvin, R. J. Schoelkopf, and L. I. Glazman, Physical Review B 86, 184514 (2012).
  * Catelani (2014) G. Catelani, Physical Review B 89, 094522 (2014).
  * Martinis et al. (2009) J. M. Martinis, M. Ansmann, and J. Aumentado, Physical Review Letters 103, 097002 (2009).
  * de Visser et al. (2014) P. J. de Visser, D. J. Goldie, P. Diener, S. Withington, J. J. A. Baselmans, and T. M. Klapwijk, Physical Review Letters 112, 047004 (2014).
  * Vool et al. (2014) U. Vool, I. M. Pop, K. Sliwa, B. Abdo, C. Wang, T. Brecht, Y. Y. Gao, S. Shankar, M. Hatridge, G. Catelani, et al., Physical Review Letters 113, 247001 (2014).
  * Wang et al. (2014) C. Wang, Y. Y. Gao, I. M. Pop, U. Vool, C. Axline, T. Brecht, R. W. Heeres, L. Frunzio, M. H. Devoret, G. Catelani, et al., Nature Communications 5, 5836 (2014).
  * Serniak et al. (2018) K. Serniak, M. Hays, G. de Lange, S. Diamond, S. Shankar, L. Burkhart, L. Frunzio, M. Houzet, and M. Devoret, Physical Review Letters 121, 157701 (2018).
  * Vepsäläinen et al. (2020) A. P. Vepsäläinen, A. H. Karamlou, J. L. Orrell, A. S. Dogra, B. Loer, F. Vasconcelos, D. K. Kim, A. J. Melville, B. M. Niedzielski, J. L. Yoder, et al., Nature 584, 551 (2020).
  * Cardani et al. (2021) L. Cardani, F. Valenti, N. Casali, G. Catelani, T. Charpentier, M. Clemenza, I. Colantoni, A. Cruciani, G. D’Imperio, L. Gironi, et al., Nature Communications 12, 2733 (2021).
  * Pop et al. (2014) I. M. Pop, K. Geerlings, G. Catelani, R. J. Schoelkopf, L. I. Glazman, and M. H. Devoret, Nature 508, 369 (2014).
  * Patel et al. (2023) H. Patel, V. Pathak, O. Can, A. C. Potter, and M. Franz, arxiv2308.02547 (2023).
  * Vool et al. (2018) U. Vool, A. Kou, W. C. Smith, N. E. Frattini, K. Serniak, P. Reinhold, I. M. Pop, S. Shankar, L. Frunzio, S. M. Girvin, et al., Physical Review Applied 9, 054046 (2018).
  * Cottet et al. (2021) N. Cottet, H. Xiong, L. B. Nguyen, Y.-H. Lin, and V. E. Manucharyan, Nature Communications 12, 6383 (2021).
  * Urwin and Arscott (1970) K. M. Urwin and F. M. Arscott, Proceedings of the Royal Society of Edinburgh Section A: Mathematics 69, 28–44 (1970).

Supplementary Material for
“A superconducting qubit based on cuprate twisted van der Waals
heterostructure”

##  Appendix A Flowermon’s spectrum

Figure 5: (a)Doublet structure of the flowermon’s spectrum as a function of
the twisting angle $\theta$ for $E_{J}=2000E_{C}$ and $E_{\kappa}=0.1E_{J}$
and $n_{g}=0$. (b) Log plot showing the behavior of the lowest doublet energy
splitting $\omega_{01}$. The red line shows the dependence of $\omega_{01}$ on
$E_{\kappa}/E_{C}$ for a fixed value of $E_{J\theta}=100\,E_{C}$ while the
blue dashed and dot-dashed lines show the dependence of $\omega_{01}$, on the
ratio $E_{J\theta}/E_{C}$ for fixed values of
$E_{\kappa}/E_{J\theta}=0.75,\,2$.

.

The flowermon’s Hamiltonian, given by Eq.(1) of main text, can be recast as
follows in the charge basis:

$H=4E_{C}\sum_{n}(n-n_{g})^{2}|n\rangle\langle
n|-\sum_{n}\left(E_{J\theta}|n\rangle\langle n+1|-E_{\kappa}|n\rangle\langle
n+2|+{\rm H.c.}\right)$ (10)

where $|n\rangle$ indicates a state of $n$ Cooper pairs. Starting from the
above Hamiltonian we calculate numerically the flowermon spectrum as a
function of the different parameters governing the Hamiltonian, namely, the
twisting angle $\theta$, the gate charge $n_{g}$, the single- and double-
Cooper pair tunneling amplitudes, $E_{J\theta}=E_{J}\cos\theta$ and
$E_{\kappa}$, and the charging energy $E_{C}$. When $\theta=45^{\circ}$, the
states with odd and even numbers of Cooper pairs are completely decoupled and
the spectrum reduces to a set of degenerate doublets, as shown in Fig.5(a).
For the low energy levels the degeneracy can be also related to the existence
of degenerate classical minima at $\phi=\pm\phi_{0}=\pm\arccos
E_{J\theta}/(4E_{\kappa})$. Indeed as hinted at by Fig.5(b) the splitting the
lowest levels scales exponentially as a function of ratio $\Delta U_{J}/E_{C}$
where $\Delta U_{J}$ is the potential barrier height defined by Eq. (2) of
main text.

In the phase basis the time-independent Schrödinger equation has the form of a
Whittaker-Hill equationUrwin and Arscott (1970)

$\frac{\partial^{2}\Psi_{n}}{\partial\varphi^{2}}+\left[a\cos\varphi-b\cos
2\varphi\right]\Psi_{n}=\lambda_{n}\Psi_{n}.$ (11)

where we introduced the normalized parameters $a=E_{J\theta}/4E_{C}$ and
$b=E_{\kappa}/4E_{C}$ and the normalized eigenvalue
$\lambda_{n}=E_{n}/4E_{C}$. Analogously to the transmon, thanks to the large
shunt capacitor, the flowermon is insensitive to gate charge fluctuations.
Specifically, it has a spectrum that becomes flat as a function of $n_{g}$
both in the limit $E_{J}\gg E_{\kappa}\gg E_{C}$ (transmon) and in the limit
$E_{\kappa}>E_{J}\gg E_{C}$ (flowermon) as shown in Fig. 6. In the two limits,
we find an exponential decrease of the charge dispersion as we decrease
$E_{C}$.

Figure 6: Structure of the low-energy spectrum of the flowermon as a function
of $n_{g}$ for different values of $E_{\kappa}/E_{C}$ and different twisting
angles.

.

##  Appendix B Quasiparticle induced decoherence

Following the route outlined by Catelani et al.Catelani et al. (2012),
starting from the microscopic Hamiltonian introduced in Ref. [Haenel et al.,
2022; Can et al., 2021] and Eqs.(5-7) of main text, we cast the quasiparticle
tunneling Hamiltonian as

$H_{T}=\sum_{{\bf kp}\sigma}t_{\bf kp}\\!\left[\left(u_{{\bf k}L}u_{{\bf
p}R}-v_{{\bf k}L}v_{{\bf
p}R}\right)\cos\frac{\varphi}{2}\gamma^{\dagger}_{{\bf k}\sigma L}\gamma_{{\bf
p}\sigma R}+i\left(u_{{\bf k}L}u_{{\bf p}R}+v_{{\bf k}L}v_{{\bf
p}R}\right)\sin\frac{\varphi}{2}\gamma^{\dagger}_{{\bf k}\sigma L}\gamma_{{\bf
p}\sigma R}{\rm+H.c.}\right]$ (12)

where the $\gamma$’s are creation and annihilation operators of Bogolyubov
quasiparticles in the left and right lead, while $u_{\bf kM}$ and $v_{\bf kM}$
denote the Bogolyubov amplitudes of lead $M$

$u_{{\bf k}M}=\sqrt{1+\frac{\xi^{2}_{\bf k}}{E^{2}_{{\bf k}M}}},\,\,v_{{\bf
k}M}=\sqrt{1-\frac{\xi^{2}_{\bf k}}{E^{2}_{{\bf k}M}}}$ (13)

where we assume $\xi_{\bf k}=\hbar^{2}k^{2}/2m-\mu$ in both leads and we
denote as usual as $E_{{\bf k}M}$ the quasiparticle energies $E_{{\bf
k}M}=\sqrt{\xi^{2}_{\bf k}+\Delta^{2}_{{\bf k}M}}$ in lead $M$ with
$\Delta_{{\bf k}M}$ indicating the corresponding gap-functions. The latter,
denoting as $\theta$ the twisting angle, can be expressed as

$\Delta_{{\bf k}L}=\Delta\cos(2\theta_{\bf k})\quad\Delta_{{\bf
k}R}=\Delta\cos(2\theta_{\bf k}-2\theta).$ (14)

Projecting the Hamiltonian (12) onto the qubit eigenstates basis and switching
to the interaction picture for the quasi-particles we get:

$\hat{V}(t)=s_{x}N_{x}(t)\sigma_{x}+c_{z}N_{z}(t)\sigma_{z}$ (15)

where

$\displaystyle\\!\\!\\!\\!\\!\\!N_{x}(t)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\sum_{{\bf kp}\sigma}\left[A^{+}_{{\bf kp}}ie^{i(E_{{\bf
k}L}-E_{{\bf p}R})t}\tilde{\gamma}^{\dagger}_{{\bf k}\sigma
L}\tilde{\gamma}_{{\bf p}\sigma R}{\rm+H.c.}\right]\\!,$ (16)
$\displaystyle\\!\\!\\!\\!\\!\\!N_{z}(t)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\sum_{{\bf kp}\sigma}\left[A^{-}_{{\bf kp}}ie^{i(E_{{\bf
k}L}-E_{{\bf p}R})t}\tilde{\gamma}^{\dagger}_{{\bf k}\sigma
L}\tilde{\gamma}_{{\bf p}\sigma R}{\rm+H.c.}\right]\\!,$ (17)
$\displaystyle{\rm and}$ $\displaystyle
s_{x}=\textrm{Tr}\left[\hat{S}\sigma_{x}\right]\,{\rm
and}\,c_{z}=\textrm{Tr}\left[\hat{C}\sigma_{z}\right]$ (18)

with the matrices $\hat{S}$ and $\hat{C}$ defined as
$\hat{S}_{ab}=\langle\psi_{a}|\sin\phi/2|\psi_{b}\rangle$ and
$\hat{C}_{ab}=\langle\psi_{a}|\cos\phi/2|\psi_{b}\rangle$ and

$A^{\pm}_{{\bf kp}}=t_{\bf kp}\left(u_{{\bf k}L}u_{{\bf p}R}\pm v_{{\bf
k}L}v_{{\bf p}R}\right).$ (19)

Note that in deriving Eq.(15) we used the approximate
$\varphi\rightarrow-\varphi$ symmetry of the flowermon eigenstates, the
symmetry becomes exact for half-integer values of $n_{g}$.

Starting from Eq. (15) we can calculate the quasiparticles-induced qubit’s
relaxation rate by means of Bloch-Redfield theory. By doing so we obtain:

$\Gamma_{\downarrow,{\rm
qp}}=\frac{s_{x}^{2}}{\hbar^{2}}\int_{-\infty}^{\infty}\\!\\!\\!\langle
N_{x}(t)N_{x}(0)\rangle e^{-i\omega_{01}t}dt$ (20)

that using Eqs.(16) and (19) can be recast as

$\displaystyle\Gamma_{\downarrow,{\rm qp}}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\frac{s_{x}^{2}}{\hbar^{2}}\sum_{{\bf
kp}}\\!\int_{-\infty}^{\infty}\\!\\!t^{2}_{\bf kp}\left(1+\frac{|\Delta_{{\bf
k}L}\Delta_{{\bf p}R}|}{E_{{\bf k}L}E_{{\bf p}R}}\right)f_{{\bf k}L}(1-f_{{\bf
p}R})\cdot$ (21) $\displaystyle\cdot e^{i(E_{{\bf k}L}-E_{{\bf
p}R}-\omega_{01})t}dt.$

Following Ref.[Yokoyama et al., 2007; Can et al., 2021] we consider the case
of coherent tunneling with in-plane momentum conservation, $t_{\bf
kp}=t\,\delta_{\bf k_{\parallel}p_{\parallel}}$, suitable to describe weakly
disordered $c$-axis junctions. This assumptions can be justified as discussed
in Ref. Bruder et al. (1995) considering the spatial dependence of the
tunneling matrix element across the barrier. A simple calculations shows that,
starting from Eq. (21) and introducing the quasiparticle density of states,
$\Gamma_{\downarrow,{\rm qp}}$ can be expressed as

$\Gamma^{\rm\theta}_{\downarrow,{\rm
qp}}=\frac{4t^{2}}{\hbar}|\langle\psi_{0}|\sin(\varphi/2)|\psi_{1}\rangle|^{2}S^{\rm\theta}_{\rm
qp}(\omega_{01})$ (22)

where we indicated as $S^{\rm\theta}(\omega_{01})$ the quasiparticle spectrum
for coherent tunneling,

Figure 7: Dependence of the first two doublets splittings on the twisting
angles. We set $E_{J}=1$THz, $E_{c}=500$MHz,$E_{\kappa}=0.1E_{J}$. Dependence
of the quasiparticle spectrum on the twisting angle. The dot indicate the
result of a numerical integration of Eq. (25) while the red dahsed line
represent the fit to an exponential
$e^{\frac{-C\Delta\sin(2\theta)^{2}}{k_{B}T}}$ with $C$ fitting coefficient.

.

$S^{\theta}_{\downarrow,{\rm qp}}(\omega_{01})=N_{0}^{2}\int\frac{d\theta_{\bf
k}}{2\pi}\int_{\Delta_{{\bf k}L}}\\!\\!\\!\\!dE\int_{\Delta_{{\bf
k}R}}\\!\\!\\!\\!dE^{\prime}\frac{EE^{\prime}+|\Delta_{{\bf k}L}\Delta_{{\bf
k}R}|}{\sqrt{E^{2}-\Delta^{2}_{{\bf
k}L}}\sqrt{{E^{\prime}}^{2}-\Delta^{2}_{{\bf
k}R}}}f(E)(1-f(E^{\prime}))\delta(E-E^{\prime}+\omega_{01}).$ (23)

where $N_{0}$ denotes the density of states at the Fermi level. It is useful
to compare the above results with case of fully incoherent tunneling, $t_{\bf
kp}=t$, applicable in the presence of strong disorder. In this case we can
write the relaxation rate as:

$\Gamma^{\rm in}_{\downarrow,{\rm
qp}}=\frac{4t^{2}}{\hbar}|\langle\psi_{0}|\sin(\varphi/2)|\psi_{1}\rangle|^{2}S^{\rm
in}_{\rm qp}(\omega_{01})$ (24)

where we indicated as $S^{\rm in}_{\rm qp}(\omega_{01})$ the quasiparticle
spectrum in the case of incoherent tunneling

$S^{\rm in}_{\rm qp}(\omega_{01})=N_{0}^{2}\int\frac{d\theta_{\bf
k}d\theta_{\bf p}}{4\pi^{2}}\int_{\Delta_{{\bf
k}L}}\\!\\!\\!\\!dE\int_{\Delta_{{\bf
p}R}}\\!\\!\\!\\!dE^{\prime}\frac{EE^{\prime}+|\Delta_{{\bf k}L}\Delta_{{\bf
p}R}|}{\sqrt{E^{2}-\Delta^{2}_{{\bf
k}L}}\sqrt{{E^{\prime}}^{2}-\Delta^{2}_{{\bf
p}R}}}f(E)(1-f(E^{\prime}))\delta(E-E^{\prime}+\omega_{01}).$ (25)

As expected the above equation yields a quasi-particle spectrum independent of
the twisting angle. Furthermore, for a $s$-wave gap, i.e. setting
$\Delta_{{\bf k}L}=\Delta_{s}$, it reduces to the result of Catelani et
al.Catelani et al. (2011) i.e.

$\Gamma^{\rm in}_{\downarrow,{\rm
qp}}\rightarrow\frac{16E_{J0}^{s}}{\hbar\,\pi}s_{x}^{2}e^{-\Delta_{s}/(k_{B}T)}e^{-\omega/(2k_{B}T)}K_{0}\left(\frac{\omega}{2k_{B}T}\right)$
(26)

where $K_{0}$ is the modified Bessel function of the second kind and we used
$E_{J0}^{s}=N_{0}^{2}t^{2}\Delta_{s}/4\pi$. A more detailed microscopic
analysis of the coherent and incoherent tunneling models and their relation
with the properties of disorder may be found in Ref.[Haenel et al., 2022].
Starting from Eq.(23), we can calculate the quasiparticle-induced relaxation
rate as a function of the quasiparticle gap $\Delta$ and the twisting angle.
By doing so we obtain the results shown in Fig. 4 of main text. Furthermore,
by fitting the numerical data, we were able to extrapolate a simple analytical
behavior for the quasiparticle spectrum, $\Gamma^{\rm\theta}_{\rm
qp}(\omega_{01})\ \sim e^{-\Delta_{d}/(k_{B}T)\sin(2\theta)^{2}}$, reminiscent
of the Yokoyama et al. result of Ref. Yokoyama et al. (2007) as shown in Fig.
7.

##  Appendix C Design and qubit parameters

The calculations in the manuscript were done using the energy ratio
$E_{J}/E_{c}=2000$. In this section, we will use realistic device parameters
to obtain this regime, remembering that $E_{c}={e^{2}}/{(2C)}$ and
$E_{J}=\Phi_{0}I_{c}/2\pi$ where $C$ is the total capacitance and $I_{c}$ is
the critical current. As a benchmark for operation at typical GHz frequencies,
we choose a Josephson energy $E_{J}$ of 1 THz and so the corresponding
critical current is $I_{c}=2$ $\mu A$. In particular, a BSCCO junction has a
critical current density $J_{c}$ of 0.1 $kA/cm^{2}$ and then the area of the
junction, i.e. the area where the two flakes are overlapped, should be of
2$\mu m^{2}$. BSCCO junctions with this area are not impossible to realize but
they are challenging because the stacking process of the two flakes must be
quick otherwise the interface will lose its superconducting properties,
becoming an insulator. Typical current junction sizes are $\approx 1000\,\mu
m^{2}$ and reducing the junction size is expected to be the most significant
limitation to the implementation of the flowermon qubit.

To achieve the corresponding capacitive energy $E_{c}=500$ MHz, so that
$E_{J}/E_{c}=2000$, we need a total capacitance $C$ of about $40fF$. This is
well within achievable parameters, and using the simulation tool Ansys we
designed and simulated a large parallel capacitor, adapting its shape to
obtain the desired value of $C$. The resulting design is displayed in Figure
8a. The capacitor thickness ($t_{c}$) is 50 nm, the distance between the two
plates of the capacitor ($g_{c}$) is 40 $\mu m$ and the total length of the
central plate ($S_{c}$) is 180 $\mu m$.

Figure 8: Control circuit. a) Parallel capacitor of $40fF$, the geometrical
dimensions of the capacitor are displayed in the figure. The capacitor
thickness ($t_{c}$) is 50 nm, the distance between the two plates of the
capacitor ($g_{c}$) is 40 $\mu m$ and the total length of the central plate
($S_{c}$) is 180 $\mu m$. b) The capacitor is coupled to a CPW $\lambda/2$
resonator with the following geometrical parameters: width of 12 $\mu m$, a
gap of 1 $\mu m$ and a length of 7 mm so that its resonance frequency is
around 8.4 GHz.The signal is transmitted to the resonator through a
transmission line (Feed line) matched to 50 $\Omega$ with a width of 5 $\mu m$
and a gap of 3 $\mu m$.

This parallel plate capacitor will be capacitively coupled with a Coplanar
Waveguide (CPW) resonator to perform the readout and to a radio-frequency line
to control the state of the qubit (Figure 8b).


# Scaling regimes for wormlike chains confined to cylindrical surfaces under
tension

Greg Morrison1,2 and D. Thirumalai3,4
1: Department of Physics, University of Houston, Houston TX 77204
2: Center for Theoretical Biological Physics, Rice University, Houston TX
77005
3: Departments of Chemistry,The University of Texas at Austin, Austin, TX,
78712
4: Department of Physics,The University of Texas at Austin, Austin, TX, 78712

###### Abstract

We compute the free energy of confinement ${\cal{F}}$ for a wormlike chain
(WLC), with persistence length $l_{p}$, that is confined to the surface of a
cylinder of radius $R$ under an external tension $f$ using a mean field
variational approach. For long chains, we analytically determine the behavior
of the chain in a variety of regimes, which are demarcated by the interplay of
$l_{p}$, the Odijk deflection length ($l_{d}=(R^{2}l_{p})^{1/3}$), and the
Pincus length ($l_{f}={k_{B}T}/{f}$, with $k_{B}T$ being the thermal energy).
The theory accurately reproduces the Odijk scaling for strongly confined
chains at $f=0$, with ${\cal{F}}\sim Ll_{p}^{-1/3}R^{-2/3}$. For moderate
values of $f$, the Odijk scaling is discernible only when ${l_{p}}\gg R$ for
strongly confined chains. Confinement does not significantly alter the scaling
of the mean extension for sufficiently high tension. The theory is used to
estimate unwrapping forces for DNA from nucleosomes.

Figure 1: Schematic diagram of the chain backbone, with a fixed distance
between monomers $a$ and resistance to bending characterized by the
persistence length $l_{p}$, the confinement to the surface of a cylinder of
radius $R$, and an external tension $f$ acting on the endpoints of the chain.

##  1 Introduction

There are compelling reasons for understanding the static and dynamics of
confined polymers because of their relevance in filtration, gel permeation
chromatography, translocation of polymers and polypeptide chains through
microporous membranes, and passage of newly synthesized proteins through the
ribosome tunnel. These and other considerations prompted several theoretical
studies, starting with pioneering studies [1, 2], which triggered subsequent
theories and simulations that probed the fate of flexible polymers in pores
with regular geometries [3, 4, 5, 6, 7, 8] as well as in the related case of
random media [9, 10] is well understood.

The situation is somewhat more complicated when considering semi-flexible
polymers or worm-like chains, WLC, in confined spaces. Spatial confinement of
WLC plays an important role in many biological systems, such as, [11, 12, 13,
14, 15], including histone wrapping in chromatin [16, 17, 18, 19] and
nanolithography [13, 20]. Here, we consider a WLC that is wrapped around a
cylinder whose radius is $R$, and is subject to mechanical force (Fig. 1).
Besides the total contour length ($L$), there are three length scales that
control the statistics of the WLC polymer. The first is the persistence length
$l_{p}$, which is a measure of the resistance of polymers to bending. For long
unconfined chains the free energy is a function of $l_{p}$. In the mean field
theory, proposed here, the polymer is globally restricted to be wound around
the cylinder, which is equivalent to restraint enforced by a soft harmonic
potential. Consequently, the Odijk length or the deflection length,
$l_{d}=(R^{2}l_{p})^{1/3}$, emerges [21] coupling the chain stiffness and
radius of confinement. In many biological contexts, the system is often under
an external field elongating the chain, such as external tension ($f$)
unravelling histone-wrapped DNA for replication [19]. An external tension or
mechanical force is captured by the Pincus length [22], $l_{f}=(\beta f)^{-1}$
with $\beta=1/k_{B}T$. The interplay of $l_{f}$, and $l_{p}$, and $l_{d}$ on
the conformations of the WLC is not fully understood. The problem has some
relevance to nucleosomes, consisting of DNA wrapped around histone proteins,
which is a building block for chromosomes. Hence, the approximate theory
developed here might illustrate an aspect of polymer theory in describing the
physics of chromatin [11, 17].

In this paper, we propose a mean-field approach to study the properties of a
wormlike chain confined to the surface of a cylinder [23] under the
application of an external tension. Because $l_{p}$ is comparable to the
contour length, $L$, excluded volume interactions may be neglected. We recover
the known [21, 24, 20, 25] dependence of free energy on the deflection length
$l_{d}=(R^{2}l_{p})^{1/3}$ at $f=0$. The theory predicts the coefficient of
the leading term of the scaling of free energy, which is in good agreement
with numerical results [23]. For moderate values of $f\neq 0$ (possibly
relevant to DNA unwrapping from nucleosomes) and strong confinement
(${l_{p}}/{R}\gg 1$), we show that that the free energy scales quadratically
with $l_{f}$. At high external tensions, we find that the effect of
confinement is perturbative, with $1-\langle
Z\rangle/L\sim\sqrt{{l_{f}}/{l_{p}}}$, independent of the radius of
confinement.

##  2 Mean Field Theory

###  2.1 Formulation and approximations

In order to understand the equilibrium properties of a cylindrically confined
WLC under tension, we developed a mean-field theory to arrive at analytically
tractable results. The energy of a configuration of the chain is determined by
a variety of terms that are sketched in Fig. 1. The system is characterized by
the spacing between monomers ($a$), the number of bonds in the chain (the
chain length $L=(N-1)a\approx Na$ with $N\gg 1$), the persistence length of
the unconfined chain ($l_{p}$), the confinement radius ($R$) aligned with the
$z$-axis, and the external tension ($f$) that is also applied in the $z$
direction. Let the position of each monomer be
${\mathbf{r}}_{i}=(x_{i},y_{i},z_{i})$, and define
${\mathbf{u}}_{i}={\mathbf{r}}_{i+1}-{\mathbf{r}}_{i}$ as the bond vector, and
$\hat{\mathbf{u}}_{i}={\mathbf{u}}_{i}/|{\mathbf{u}}_{i}|$ the unit bond
vector. The statistics of a surface-confined chain can be described by a
constrained Kratky-Porod (KP) Hamiltonian[26] $\beta
H_{KP}=\frac{l_{p}}{a}\sum_{n=0}^{N}\hat{\mathbf{t}}_{i}\cdot\hat{\mathbf{t}}_{i+1}-\beta
f(z_{N}-z_{0})$, with confinement to the surface of the cylinder requiring
$x_{i}^{2}+y_{i}^{2}=R^{2}$ for all $i$.

The KP model is mathematically difficult to work with because of two rigid
constraints: the fixed bond length, $|{\mathbf{u}}_{i}|=a$, and the constraint
that monomers be spatially constrained transverse to the $z$ axis, which
should be enforced through the relation, $x_{i}^{2}+y_{i}^{2}=R^{2}$. On the
mean field level, we replace these rigid constraints for the confined WLC with
softer harmonic restraints (an approach that has been fruitfully applied
previously [27, 25, 28, 29, 30]). The form of the MF Hamiltonian can be found
by writing the monomer distribution function explicitly, in order to identify
a physically meaningful harmonic approximation to the rigid constraints. It is
straightforward to show that (up to a constant) the statistical weight is
$\Psi_{S}=(\prod_{n=0}^{N}\delta[x_{n}^{2}+y_{n}^{2}-R^{2}])\times(\prod_{n=1}^{N}\delta[{|{\Delta{\mathbf{r}}}_{n}|-a}])\times(\prod_{n=1}^{N-2}e^{-l_{p}/2a^{3}({\Delta{\mathbf{r}}}_{n+1}-{\Delta{\mathbf{r}}}_{n})^{2}})\times
e^{+\beta f(z_{L}-z_{0})}$ where
${\Delta{\mathbf{r}}}={\mathbf{r}}_{n+1}-{\mathbf{r}}_{n}$ and
$\beta=1/k_{B}T$. The first term enforces confinement of the chain to the
cylinder (affecting only the $x$ and $y$ coordinates of the polymer), the
second term enforces the constant monomer spacing constraint, the third term
accounts for the WLC’s resistance to bending, and the fourth term accounts for
the external tension. Each of the $\delta$ functions may be written as
$\delta(x_{n}^{2}+y_{n}^{2}-R^{2})\propto\int
dk_{n}e^{iak_{n}[(x_{n}^{2}+y_{n}^{2})/R^{2}-1]}$ and
$\delta(|{\Delta{\mathbf{r}}}|-a)\propto\int
d\lambda_{n}e^{-ia\lambda_{n}[{\Delta{\mathbf{r}}}^{2}_{n}/a^{2}-1]}$, leading
to,

$\displaystyle\Psi_{S}$ $\displaystyle\propto$
$\displaystyle\int_{-i\infty}^{i\infty}\prod_{n}dk_{n}d\lambda_{n}\exp\bigg{[}-a\sum_{n=1}^{N}k_{n}\bigg{(}\frac{x_{n}^{2}+y_{n}^{2}}{R^{2}}-1\bigg{)}$
$\displaystyle-a\sum_{n=1}^{N-1}\lambda_{n}\bigg{(}\frac{{\Delta{\mathbf{r}}}_{n}^{2}}{a^{2}}-1\bigg{)}-\beta
f(z_{L}-z_{0})-\frac{1}{2}a\sum_{n=1}^{N-2}l_{p}\frac{({\Delta{\mathbf{r}}}_{n+1}-{\Delta{\mathbf{r}}}_{n})^{2}}{a^{4}}\bigg{]}.$

The $\delta$ functions constrain enforce a fixed monomer separation and fixed
transverse distance in the confined dimension. Here, we focus solely on
surface confinement, but note that a MF approach has been applied to the
confinement to the interior and surface [27] of a sphere. An equivalent volume
constraint could be applied the interior of the cylinder through a constraint
of $x_{i}^{2}+y_{i}^{2}\leq R^{2}$. However, a harmonic constraint on volume
confinement requires an estimate of the average distance of each monomer in
the transverse direction (which cannot be predicted at the mean field level).

Although exact, the formulation in Eq. 2.1 is difficult to work with directly.
Analytical progress becomes possible by assuming that the integrals over the
Fourier variables are sharply peaked, in the same manner as in our previous
studies [27, 31]. The partition function
$Z=\prod_{n}d^{3}{\mathbf{r}}_{n}\Psi_{S}(\\{{\mathbf{r}}_{n}\\})\equiv\int\prod_{n}d\lambda_{n}dk_{n}\exp(-{\cal{F}}[\\{\lambda_{n},k_{n}\\}])$
defining the nondimensional free energy functional ${{\cal{F}}}$ as an
integral over all the monomer coordinates, and the linearity of the Fourier
transform allows us to write
${\cal{F}}={\cal{F}}_{x}+{\cal{F}}_{y}+{\cal{F}}_{z}-a\sum_{n=1}^{N}k_{n}-a\sum_{n=1}^{N-1}\lambda_{n}$,
with
$e^{-{\cal{F}}_{x}}\equiv\int\prod_{n}dx_{n}e^{-{\cal{H}}_{c}[\\{x_{n}\\}]}$
and
$e^{-{\cal{F}}_{z}}\equiv\int\prod_{n}dz_{n}e^{-{\cal{H}}_{u}[\\{z_{n}\\}]}$,
where we define the confined Hamiltonian as, ${\cal{H}}_{c}[\\{x_{n}\\}]\equiv
a\sum_{n=1}^{N-2}{l_{p}(\Delta x_{n+1}-\Delta
x_{n})^{2}}/{2a^{4}}+a\sum_{n=1}^{N-1}\lambda_{n}{\Delta
x^{2}_{n}}/{a^{2}}+a\sum_{n=1}^{N}k_{n}{x_{n}^{2}}/{R^{2}}$. Similarly, the
unconfined Hamiltonian is, ${\cal{H}}_{u}[\\{z_{n}\\}]\equiv
a\sum_{n=1}^{N-2}{l_{p}(\Delta z_{n+1}-\Delta
z_{n})^{2}}/{2a^{4}}+a\sum_{n=1}^{N-1}\lambda_{n}{\Delta
z^{2}_{n}}/{a^{2}}+\beta f\sum_{n=1}^{N-1}\Delta z_{n}.$ If the free energy is
sharply peaked around some
$\\{\lambda_{n},k_{n}\\}=\\{\lambda_{n}^{*},k_{n}^{*}\\}$ that minimizes
${\cal{F}}$, the partition function can be written approximately as $Z\approx
Z^{*}=e^{-{\cal{F}}^{*}}$. Because the Hamiltonians ${\cal{H}}_{c}$ and
${\cal{H}}_{u}$ are uncoupled and quadratic in the monomer coordinates, it is
straightforward to integrate over the internal coordinates of the polymer
exactly. Along the confined axes, we find
${\cal{F}}_{x}={\cal{F}}_{y}=\frac{1}{2}\log[{\mbox{Det}}({\mathbf{Q}})]$,
where the elements of the matrix $({\mathbf{Q}})_{nm}$ are the coefficients
associated with $x_{n}x_{m}$ in Eq. 2.1. The explicit form of
$({\mathbf{Q}})_{nm}$ is given in Appendix A of an earlier work [27]. In the
unconfined direction, completing the square and noting the translational
invariance of the system along the cylinder axis, allows us to write
${\cal{F}}_{z}=\frac{1}{2}\log[{\mbox{Det}}({\mathbf{P}})]e^{\frac{1}{4}(\beta
f)^{2}a\sum_{n=1}^{N-1}\lambda_{n}^{-1}}$, with [31]
$({\mathbf{P}})_{ij}=\lambda_{i}\delta_{ij}-l_{p}\delta_{i,j\pm 1}/2a^{2}$.
These expressions can in principle be used to determine the stationary phase
values for the Fourier variables by setting
$\partial{\cal{F}}/\partial\lambda_{n}|_{\\{\lambda_{n},k_{n}\\}=\\{\lambda_{n}^{*},k_{n}^{*}\\}}$
= $\partial{\cal{F}}/\partial
k_{n}|_{\\{\lambda_{n},k_{n}\\}=\\{\lambda_{n}^{*},k_{n}^{*}\\}}=0$. However,
the resulting $2N-1$ equations become intractable for large $N$, thus making
it necessary to make additional approximations.

The matrices ${\mathbf{P}}$ and ${\mathbf{Q}}$ are bidiagonal and tridiagonal,
respectively, with a regular structure except near the endpoints of the chain.
The high symmetry of the matrices that underly the equations suggest that we
should seek symmetric solutions for the stationary values of $\lambda_{n}^{*}$
and $k_{n}^{*}$. In both the unconfined[31] and spherically confined[27]
cases, it was shown that the tractable equations that reproduced exact
theoretical results could be found for mean field parameters separated into
bulk terms and endpoint terms[31, 27], with $\lambda_{n}^{*}=\lambda$ and
$k_{n}=k$ for the interior points on the chain (those with $2<n<N-2$). Excess
endpoint fluctuations generally require that $\lambda_{1}=\lambda_{N-1}\neq 1$
and $k_{1}=k_{N}\neq k\neq k_{2}=k_{N-1}$ (with the equalities due to symmetry
arguments). This approximation allows us to write
${\cal{H}}_{c}[\\{x_{n}\\}]={\cal{H}}_{c}^{(b)}[\\{x_{n}\\}]+{\cal{H}}_{c}^{(e)}[\\{x_{n}\\}]$
and ${\cal{H}}_{u}[\\{\Delta z_{n}\\}]={\cal{H}}_{c}^{(b)}[\\{\Delta
z_{n}\\}]+{\cal{H}}_{c}^{(e)}[\\{\Delta z_{n}\\}]$, where the superscript $b$
denotes an extensive bulk term, depending on the mean field parameters
$\lambda$ and $k$. The superscript $e$ denoting an intensive endpoint term
depending on the more complicated values of the mean field variables near the
ends of the chain.

In the continuum limit, it can be shown that the bulk form of the Hamiltonians
become,

$\displaystyle{\cal{H}}_{c}^{(b)}[x(s)]$ $\displaystyle=$
$\displaystyle\int_{0}^{L}ds\bigg{(}\frac{l_{p}}{2}\ddot{x}^{2}(s)+\lambda\dot{x}^{2}(s)+k\frac{x^{2}(s)}{R^{2}}\bigg{)}$
(2) $\displaystyle{\cal{H}}_{u}^{(b)}[\dot{z}(s)]$ $\displaystyle=$
$\displaystyle\int_{0}^{L}ds\bigg{(}\frac{l_{p}}{2}\ddot{z}^{2}(s)+\lambda\dot{z}^{2}(s)-\beta
f\dot{z}(s)\bigg{)}.$ (3)

The form for ${\cal{H}}_{c}^{(b)}$ is identical to the spherically confined
Hamiltonian of the wormlike chain in our previous work [27], while the form of
${\cal{H}}_{u}^{(b)}$ is identical to that for an unconfined chain [31]. Note
that $\lambda$ is the same in the two Hamiltonians in Eq. 3 because of the
condition $\langle{\mathbf{u}}^{2}\rangle=\langle
u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\rangle=1$ couples the three components of the
bending vector. However, the mean field parameter $k$ occurs only in the
confined Hamiltonian and enforces $\langle x^{2}+y^{2}\rangle=R^{2}$
constraint. The assumption that $\lambda$ is isotropic causes inaccuracies in
the predicted mean extension of a chain using the MF approach [32]. A MF
theory that avoids this assumption has been developed [30], which does at the
cost of greater complexity in the model. As we will show below, the expected
scaling of the free energy in the limits of high tension and strong
confinement are both recovered using an isotropic assumption, suggesting that
the overall scaling laws will be accurate but with potentially inaccurate
coefficients. We also ignore the endpoint effects by neglecting the
Hamiltonians ${\cal{H}}_{c}^{(e)}$ and ${\cal{H}}_{u}^{(e)}$. This
approximation simplifies the mathematics greatly, but restricts our analysis
to very long chains. Neglecting the endpoint effects, the mean field equations
for the extensive contribution to the free energy in the continuum are
$\partial{\cal{F}}/\partial\lambda=\partial{\cal{F}}/\partial k=0$, with
$e^{-{\cal{F}}}=\int{{{\cal{D}}}}[{\mathbf{r}}(s)]e^{-{\cal{H}}_{c}^{(b)}[x(s)]-{\cal{H}}_{c}^{(b)}[y(s)]-{\cal{H}}_{u}^{(b)}[z(s)]+L(\lambda+k)}$.

###  2.2 Calculation of the Free Energy

It is straightforward to perform the path integrals over the confined [27] and
unconfined [31] dimensions to explicitly compute the partition function
$e^{-{\cal{F}}}$, from which it is possible to derive the mean field equations
for constant $\lambda$ and $k$. One can readily recognize that Eq. 3 describes
a quantum harmonic oscillator after a change of variables, for which an exact
propagator is known[33]. In the confined dimension, it is straightforward to
show[27] that the action in Eq. 2 is minimized by a path satisfying
$\frac{l_{p}}{2}x^{(4)}(s)-\lambda\ddot{x}(s)+\frac{k}{R^{2}}x(s)=0$. It is
readily observed (after some tedious mathematics) that the solution is
expressible in terms of the frequencies
$w_{\pm}^{2}=\frac{\lambda}{l_{p}}(1\pm\sqrt{1-{2kl_{p}}/{\lambda^{2}R^{2}}}\
)$. It is possible to integrate over the internal degrees of freedom of the
chain exactly [27], which results in an unwieldy expression. However, with the
assumption that $\sinh(L\omega_{\pm})\approx\cosh(L\omega_{\pm})\approx
e^{L\omega_{\pm}}/2$ (satisfied for $L/l_{p}\gg 1$), it is possible to
simplify the confinement free energy as,

$\displaystyle\frac{{\cal{F}}}{L}\approx\bigg{(}\omega_{+}+\omega_{-}-k\bigg{)}+\bigg{(}\sqrt{\frac{\lambda}{2l_{p}}}-\lambda-\frac{(\beta
f)^{2}}{4\lambda}\bigg{)},$ (4)

where the first term is the contribution from integration over the chain
configurations along the two confined axes, and the second is the contribution
from the single unconfined axis. Endpoint effects are neglected in deriving
Eq. 4, and is only valid for very long chains (where $L$ is larger than all
other length scales in the problem). The average extension of the chain, under
tension, can be calculated using,

$\displaystyle\langle Z\rangle=\frac{\partial{\cal{F}}}{\partial(\beta
f)}=\frac{\beta fL}{2\lambda(l_{p},R,f)},$ (5)

where we have explicitly included the dependence of the mean field solution of
$\lambda$ on the physical parameters at play. Eq. 5 is similar to the result
found in the case of an unconfined chain under an external tension[28], which
is straightforward to evaluate once the mean field solution for $\lambda$ is
known.

To determine the free energy and mean extension, we must solve the MF
equations $\partial{\cal{F}}/\partial k=\partial{\cal{F}}/\partial\lambda=0$.
The resulting equations are greatly simplified by noting that

$\displaystyle k=\frac{l_{p}R^{2}}{2}\omega_{1}^{2}\omega_{2}^{2}\qquad$
$\displaystyle\qquad\lambda=\frac{l_{p}}{2}(\omega_{1}^{2}+\omega_{2}^{2})$
(6)
$\displaystyle\frac{\partial\omega_{\pm}}{\partial\lambda}=\pm\frac{\omega_{\pm}}{l_{p}(\omega_{+}^{2}-w_{-}^{2})}\quad$
$\displaystyle\quad\frac{\partial\omega_{\pm}}{\partial
k}=\mp\frac{1}{l_{p}R^{2}\omega_{\pm}(\omega_{+}^{2}-w_{-}^{2})}.$ (7)

After some algebra, the variational equations for $\lambda$ and $k$,
respectively, become

$\displaystyle\frac{1}{l_{p}(\omega_{+}+\omega_{-})}=1-\frac{1}{\sqrt{8\lambda
l_{p}}}-\frac{(\beta f)^{2}}{4\lambda^{2}}\quad$
$\displaystyle\quad\frac{1}{l_{p}(\omega_{+}+\omega_{-})}=R^{2}\omega_{+}\omega_{-}.$
(8)

As the left-hand side of both equalities in the above equation are identical,
we can readily solve for $k$ in terms of $\lambda$, with,

$\displaystyle k=\frac{l_{p}}{2R^{2}}\bigg{(}1-\frac{1}{\sqrt{8\lambda
l_{p}}}-\frac{1}{4}\left(\frac{\beta f}{\lambda}\right)^{2}\bigg{)}^{2}$ (9)

so that the confinement parameter $k$ at the mean field level can be
determined exactly. The mean field parameter $\lambda$, enforcing the
inextensibility of the chain, requires the solution of a complicated equation
in Eq. 8, after substitution of the exact solution for $k$. Although it may
not be possible to solve for the exact values analytically for all $l_{p}$,
$R$, and $f$, it is straightforward to numerically determine the mean field
values accurately for any $l_{p}$ and $R$. The asymptotic behavior of the
roots can be readily determined in certain limits, as discussed in the next
section.

##  3 Results

###  3.1 Scaling for weak confinement and weak tension

An asymptotic solution to the Mean Field equations can be determined in a
variety of parameter regimes. The simplest scenario is the limit where both
the cylindrical confinement and the external tension are weak. In the limit of
$R/l_{f}\ll 1$ (with $l_{f}=k_{B}T/f$ the Pincus length) and $l_{p}/R\ll 1$,
the chain is weakly confined and we expect the mean field solution to be a
perturbation on the solution to the unconfined three dimensional mean field
theory. It has been shown [31] that $\lambda\sim 9/8l_{p}$ in the absence of
external force or confinement. In order to determine the asymptotic behavior
of the mean field parameter, we expand in a series for small values of the
rescaled variables $\tilde{l}_{p}\equiv l_{p}/R\ll 1$ and $\varphi=R/l_{f}\ll
1$, with $\lambda
R=\frac{9}{8\tilde{l}_{p}}+\sum_{m,n=0}^{\infty}b_{nm}{\tilde{l}_{p}}{}^{n}\varphi^{m}$
for some coefficients $b_{nm}$. Substitution for the exact solution for $k$ in
Eq. 9 into the MF equation for $\lambda$ in eq. 8 allows us to perform a
series expansion in higher order of $\tilde{l}_{p}$ and $\varphi$, both
assumed to be small. Iteratively solving for the lowest order coefficients
$b_{ij}$ shows
$\lambda\approx\frac{9}{8l_{p}}-\frac{4l_{p}}{9R^{2}}+\frac{4l_{p}}{9l_{f}^{2}}+O(l_{p}^{3}/R^{4})$,
to leading order in $\tilde{l}_{p}$ and $\varphi$. To leading order, the free
energy under weak confinement and weak external tension is

$\displaystyle\frac{{\cal{F}}}{L}\approx\frac{9}{8l_{p}}+\frac{2l_{p}}{9R^{2}}\left(1-\frac{8l_{p}^{2}}{81R^{2}}\right)-\frac{2l_{p}}{3l_{f}^{2}}\qquad(l_{p}\ll
R,R/l_{f}f\ll 1).$ (10)

Note that retaining higher order terms in $\lambda$ does not affect the
scaling coefficients in eq. 10. It is interesting to note that the leading
order contribution of the tension is independent of the radius of the
cylinder, with the confinement only entering into the free energy through
coupling with the persistence length. This is due to the distinct axes over
which each energetic contribution acts, each of which are perturbative in this
limit. Not surprisingly for weakly confined chains, the deflection length
$l_{d}=(l_{p}R^{2})^{1/3}$ does not enter into the free energy. We note that
these scaling coefficients are not likely to be precise, as has been
previously noted for the MF solutions in multiple contexts [28, 30, 32].
However, the scaling with each variable is expected to be accurate.

For a weakly confined chain, the average extension in eq. 5 becomes

$\displaystyle\frac{\langle
Z\rangle}{L}\approx\frac{4}{9}\frac{l_{p}/l_{f}}{1-{2}\left(\frac{4l_{p}}{9R}\right)^{2}},$
(11)

to fourth order in $\tilde{l}_{p}$ and $\varphi$, growing linearly with $f$ as
long as $l_{p}\ll R$. The linear increase for low forces and weak confinement
is shown in the purple triangles of Fig. 2 (A), satisfying the expected linear
scaling.

Figure 2:  (A) Average chain extension (Eq. 5) as a function of $f$ (with
$R=1$ held fixed) for various values of $l_{p}$: $l_{p}/R=1$ (purple
triangles), 10 (red circles), and 100 (blue squares). The extension increases
linearly with force for small $f$, and the linear scaling is unperturbed by
confinement effects. (B) The approach to full extension, $1-\langle
Z\rangle/L$, as a function of $f$ scales as $f^{-1/2}$ for large external
tension, but deviations from this scaling occur for very stiff chains.

###  3.2 Strongly confined chains under weak tension

It is also possible to determine the scaling behavior of the polymer under
strong confinement (with $l_{p}\gg R$) but still constraining the external
tension to be weak (with $R\beta f\ll 1$). A one-dimensional WLC in the
absence of tension on the mean field level will satisfy $\lambda\sim
1/8l_{p}$, and we expect that $\lambda$ must converge to this value for
sufficiently small $R$ or sufficiently large $l_{p}$ (since transverse
fluctuations must vanish in either limit). The effect of confinement is
contained in the mean field variable $k$, which should capture the transverse
statistics of the chain. It is known that the Odijk length scale,
$l_{d}\propto(R^{2}l_{p})^{1/3}$, emerges for strongly confined WLCs to the
interior of a cylinder [21, 25]. A stiff chain will predominantly be aligned
with the cylinder axis, and $l_{d}$ would be the typical distance along the
chain for transverse fluctuations (a consequence of the mean field approach)
to encounter the walls (causing a deflection). A similar scaling has been
observed for surface confined chains [23] and we expect the deflection length
to emerge at the MF level because the hard constraint is replaced by a soft
harmonic potential. As a consequence, the stiff chain ‘deflects’ off the soft
harmonic potential rather than a rigid wall.

In the limit of $\tilde{l}_{p}\equiv l_{p}/R\gg 1$ and $\varphi=R/l_{f}\ll 1$
we expect the leading order behavior of $\lambda$ can be recovered with the
ansatz
$\lambda\approx\frac{1}{8l_{p}}+\sum_{n=4}^{\infty}c_{n}\tilde{l}_{d}^{-n/3}+\sum_{m=1}^{\infty}d_{m}\varphi^{m}$.
This is similar to the expression for $\lambda$ in the weakly confined case
but with the expansion in terms of small $\varphi$ and large
$\tilde{l}_{p}^{1/3}$, and ignores cross-terms (assumed to be higher order, as
was the case for weak confinement and low tension). The restriction of $n\geq
4$ in the sum over $\tilde{l}_{p}$ ensures we recover the expected 1-d scaling
for $R=0$, since the leading order behavior is expected to be
$\lambda=1/8l_{p}$ for $R=0$. Substitution into Eqs. 8 using 9 and taking the
limit of small $\varphi$ and large $\tilde{l}_{p}$ shows that $c_{4}=0$ and
$c_{5}=1/(4\times 2^{1/3})$ cancel the lowest order terms in $\tilde{l}_{p}$,
while $d_{1}=0$ and $d_{2}=4$ cancels the leading term in $\varphi$.
Substitution of this solution into the free energy yields

$\displaystyle\frac{{\cal{F}}}{L}$ $\displaystyle\approx$
$\displaystyle\frac{3}{2^{5/3}l_{d}}-\frac{2l_{p}}{l_{f}^{2}}+\frac{l_{d}}{l_{f}^{2}}+\frac{1}{8l_{p}},$
(12)

with the first term the leading order in the confinement, the second to
leading order in $l_{f}$, the third term the lowest order coupling between the
competing length scales $l_{d}$ and $l_{f}$, and the fourth term the
1-dimensional MF free energy. In the absence of force, this expression for the
free energy at a large eternal tension scales as $L/l_{d}$, agreeing with the
known result for a strongly confined WLC [21]. The coefficient of the leading
term in Eq. 12 is $3/2^{5/3}\approx 0.945$, which will dominate in the limit
of $R\to 0$. This compares reasonably well with the theoretical prediction for
a surface-confined chain, with Eq. 38 in a previous study [23], which predicts
a coefficient of 0.8416 by solving an exact nonlinear Fokker-Plank equation.
We also note that the coefficient of the leading term of Eq. 12 is exactly
twice the value Eq. [8] in an earlier report of DNA in confined to the
interior of a cylinder [25]. This agreement along with the accurate prediction
of the scaling coefficient in the previously studied case of spherical
confinement using the mean field theory [27] suggest that even the calculated
numerical coefficients arising may generally be reliable.

Figure 3: (A) Numerical roots for the extensive component of the mean field
free energy for $R\beta f=0$ (black stars), 1 (purple triangles), 10 (red
circles), and 100 (blue squares) as $l_{p}$ is varied and $R$ held fixed. The
scalings derived for the three asymptotic regimes are indicated by the gray
dashed lines. The free energy is shifted by a factor $L/l_{f}$ to counter the
leading term in Eq. 15. (B) The coefficients $\alpha$ for the curves in (A)
show that the scaling of ${\cal{F}}\sim l_{p}^{-1}$ and ${\cal{F}}\sim
l_{p}^{-1/3}$ are recovered for sufficiently weak or sufficiently strong
confinement respectively. Furthermore, ${\cal{F}}\sim l_{p}^{-1/2}$ is
observed only for very large forces.

The scaling is confirmed by the numerical calculation for $l_{p}/R\gtrsim 10$
at $f=0$ in the black stars in Fig. 3, and for non-zero forces the scaling is
recovered for sufficiently strong confinement. Fig. 3(B) shows that non-zero
tensions can significantly alter the onset of the Odijk regime. Even an
external tension as low as, ${R}/{l_{f}}=1$ (the purple triangles), delays the
onset of the transition to the Odijk scaling,
$\alpha=\partial\log({\cal{F}})/\partial\log(l_{p}/R)=-\frac{1}{3}$, by orders
of magnitude in $l_{p}/R$. In particular, we note that modeling the histone as
a long cylinder of radius $R\approx 3.15$nm and DNA with $l_{p}\approx 50$nm
will have a free energy of confinement satisfying the Odijk scaling at $f=0$
(since $l_{p}/R\approx 16$). The binding of DNA to histones is known to be
interrupted by forces on the order of $f\approx 10$pN [19] and the unbinding
of DNA already bound to histones occurs [19, 18] at forces $\approx 20-30$pN
(or $7.7\lesssim\frac{R}{l_{f}}\lesssim 23.0$). The intermediate external
tension in the red squares of Fig. 3 (with ${R}/{l_{f}}=10$) falls within this
range, indicating that the Odijk scaling may not be discernible by tensile
forces over a wide range of biologically relevant conditions. It is
interesting to note that confinement does not have a similarly strong affect
on the scaling for the extension of a WLC under tension. The extension of the
chain for small $f$ and large $l_{p}/R$ becomes

$\displaystyle\frac{\langle
Z\rangle}{L}\approx\frac{4l_{p}/l_{f}}{1+\left(\frac{2R}{l_{p}}\right)^{2/3}},$
(13)

with the variations in the denominator being weak due to the approximation
$l_{p}/R\gg 1$. Note that the scaling is linear in $l_{p}$, despite the
emergence of the deflection length as the dominant length scale in the free
energy in Eq. 12. We expect that even stiff chains will deviate only slightly
from linear scaling (confirmed in the red squares and blue circles of Fig.
2(B). The leading term in Eq. 13 differs from the extension of a weakly
confined chain by a factor of $\frac{1}{9}$, which implies that strongly
confined chains are expected to have a greater average extension at the same
external tension as weakly confined chains (a physically sensible result).
This is confirmed in Fig. 4, which shows that a decrease in the radius of the
cylinder at fixed $l_{p}$, decreases the midpoint of extension ($f_{m}$, the
force at which $\langle Z\rangle/L=0.5$) from $l_{p}\beta f_{m}\sim 2$ by
nearly an order of magnitude. For parameters matching the wrapping of DNA
around a histone core, with $l_{p}\approx 16$ and $R\approx 3.15$nm, we
predict that the midpoint of the transition is greatly reduced from its
unconfined value, with $f_{m}(R\to\infty)\approx 0.16$pN for the unconfined
chain[28] and $f_{m}(R=3.15)\approx 0.025$pN.

Figure 4: Extension $\langle Z\rangle/L$ for fixed $l_{p}=1$ as $R$ and $f$
are varied. The midpoint of the force extension curve for sufficiently large
$R/l_{p}\sim 10^{3}$ occurs near $l_{p}\beta f_{m}\sim 2$, but is reduced by
nearly a factor of 8.8 for small $R/l_{p}\sim 10^{-3}$ to $l_{p}\beta
f_{m}\sim 0.23$, in good agreement with the predicted reduction by a factor of
9 between weakly and strongly confined chains.

###  3.3 Confined chains under high tension

Finally, we consider the limit of external tensions that are strong enough to
be comparable to confinement effects, satisfying $\frac{R}{l_{f}}\gtrsim
l_{p}/R$. To determine the asymptotic form of the free energy in this limit,
we define $\tilde{l}_{p}=l_{p}/R\equiv\sigma\times R/l_{f}=\sigma\varphi$ for
some unconstrained $\sigma$ (assumed finite, but not restricted to be large or
small), and expand $\lambda$ in the limit of large $\varphi$. In this limit,
it is straightforward to show that $\lambda
R\approx\frac{\varphi}{2}+\frac{3}{8\sqrt{\sigma}}+\frac{(9\sigma-8)}{32\sigma\varphi}$
to leading order in $\varphi\gg 1$. In the high-force and high-tension limit,
the free energy becomes,

$\displaystyle\frac{{\cal{F}}}{L}$ $\displaystyle\approx$
$\displaystyle-\frac{\varphi}{R}+\frac{3}{2R\sqrt{\sigma}}+\frac{1}{2R\varphi}\left(1+\frac{9}{16\sigma}\right)$
(14) $\displaystyle=$ $\displaystyle-\beta f+\frac{3}{2}\sqrt{\frac{\beta
f}{l_{p}}}+\frac{9}{32l_{p}}+\frac{1}{2\beta fR^{2}}.$ (15)

Note that the free energy here does not depend on the deflection length
$l_{d}\sim(l_{p}R^{2})^{1/3}$ that was seen for strongly confined chains with
low tension. Thus, a strong force significantly alter the free energy of
confinement. The predicted scaling is confirmed in the blue circles in Fig. 3
for very large forces (with $R\beta f=100$). We do not observe the expected
limits of ${\cal{F}}\sim(l_{p}/l_{f})^{-1/2}$ or
${\cal{F}}\sim(l_{p}R^{2})^{-1/3}$ for strong confinement and intermediate
tension (the purple triangles or red squares with ${R}/{l_{f}}=1$ or 10
respectively), but instead an extended crossover region between the expected
scaling laws emerges. The average extension of the chain in this limit is
similar to the unconfined WLC, with

$\displaystyle\frac{\langle Z\rangle}{L}=1-\frac{3}{4\sqrt{l_{p}\beta f}},$
(16)

with the highest order dependence on the confinement radius scaling as
$(\frac{R}{l_{f}})^{-2}$, which is assumed to be small in this limit. Note
that the inextensibility of the chain is respected in the mean field theory
(since $\langle Z\rangle\to L$ as $f\to\infty$), and that Eq. 16 (valid for
sufficiently high forces) is identical to that of the unconfined chain. The
approach to full extension is shown in Fig. 2(B), confirming the scaling of
$1-\langle Z\rangle/L\sim f^{-1/2}$ for large $f$ regardless of the strength
of confinement. Deviations from this scaling occur only for much weaker forces
(with $l_{p}\beta f\lesssim 1$) for strongly confined chains.

##  4 Conclusions

In this paper, we calculated the free energy and linear extension of a
wormlike chain confined to the surface of a cylinder with an applied external
tension using a mean field approach. We conclude with the following additional
comments. (1) Our method recovers the Odijk scaling of the free energy in the
absence of force and the one-dimensional extension profile for a WLC in the
limit of small cylinder radius. (2) The coefficient in Eq. 12 in the free
energy expression obtained here is fairly close to the result obtained
previously [23], which shows that the mean field theory not only predicts the
scaling relation but also yields the coefficients that are fairly accurate.
(3) For a nucleosome, $\frac{l_{p}}{R}$ lies in the range (5-10) depending on
the concentration and valence of counter ions. Moreover, experiments [18] show
that the first stage of DNA unwrapping occurs at $f\approx(3-5)$ pN. For these
parameters, the theory predicts that
$\alpha=\partial\log({\cal{F}})/\partial\log(l_{p}/R)\approx-0.5$ (red curve
in Fig. 3B), which corresponds to strong external tension and modest
confinement.

Acknowledgements: DT is grateful to Fyl Pincus for free lessons, with
historical connotations, on polymer physics. This work was supported by grants
from the National Science Foundation (CHE 2320256, PHY 2014141 and PHY
2019745) and the Welch Foundation through the Collie-Welch Chair (F-0019).

## References

  * [1] E.F. Casassa and Y Tagami.  An equilibrium theory for exclusion chromatography of branched and linear polymer chains.  Macromolecules, 2:14–26, 1969.
  * [2] M Daoud and P.G. de Gennes.  Statistics of macromolecular solutions trapped in small pores.  J. de. Physique, 38:85–93, 1977.
  * [3] F. Brochard-Wyart and E. Raphael.  Scaling theory of molten polymers in small pore.  Macromolecules, 23:2276–2280, 1990.
  * [4] G. Morrison and D. Thirumalai.  The shape of a flexible polymer in a cylindrical pore.  J. Chem. Phys., 122:194907, 2005.
  * [5] S. Jun, D. Thirumalai, and BY Ha.  Compression and stretching of a self-avoiding chain in cylindrical nanopores.  Phys. Rev. Lett, 101:138101, 2008.
  * [6] Cristian Micheletti, Davide Marenduzzo, and Enzo Orlandini.  Polymers with spatial or topological constraints: Theoretical and computational results.  Physics Reports, 504(1):1–73, 2011.
  * [7] Youngkyun Jung, Chanil Jeon, Juin Kim, Hawoong Jeong, Suckjoon Jun, and Bae-Yeun Ha.  Ring polymers as model bacterial chromosomes: confinement, chain topology, single chain statistics, and how they interact.  Soft Matter, 8(7):2095–2102, 2012.
  * [8] E Werner and B Mehlig.  Confined polymers in the extended de gennes regime.  Physical review E, 90(6):062602, 2014.
  * [9] S. F. Edwards and D. Thirumalai.  The size of a polymer in random media.  J. Chem. Phys., 89:2435–2441, 1988.
  * [10] D. Thirumalai.  Isolated polymer molecule in a random environment.  Phys. Rev. A, 37:269–276, 1988.
  * [11] Assaf Amitai and David Holcman.  Polymer physics of nuclear organization and function.  Physics Reports, 678:1–83, 2017.
  * [12] Jeff ZY Chen.  Theory of wormlike polymer chains in confinement.  Progress in Polymer Science, 54:3–46, 2016.
  * [13] Walter Reisner, Jonas N Pedersen, and Robert H Austin.  DNA confinement in nanochannels: physics and biological applications.  Reports on progress in physics, 75(10):106601, 2012.
  * [14] Bae-Yeun Ha and Youngkyun Jung.  Polymers under confinement: single polymers, how they interact, and as model chromosomes.  Soft Matter, 11(12):2333–2352, 2015.
  * [15] Omar A Saleh.  Perspective: Single polymer mechanics across the force regimes.  The Journal of chemical physics, 142(19), 2015.
  * [16] J. Langowski.  Polymer chain models of dna and chromatin.  Eur. Phys. J. E, 19:241, 2006.
  * [17] R. Everaers and H. Schiessel.  The physics of chromatin.  J. Phys. Cond. Mat., 27:060301, 2015.
  * [18] Y Cui and C Bustamante.  Pulling a single chromatin fiber reveals the forces that maintain its higher-order structure.  Proc. Natl. Acad. Sci., 97:127, 2000.
  * [19] M. L. Bennik, S. H. Leuba, G. H. Leno, J. Zlatanova, B. G. de Grooth, and J. Greve.  Unfolding individual nucleosomes by stretching single chromatin fibers with optical tweezers.  Nat. Struct. Biol., 8:606–610, 2001.
  * [20] W. Reisner, K. J. Morton, R. Riehn, Y.M. Wang, Z. Yu, M. Rosen, J. C. Sturm, S. Y. Chou, E. Frey, and R. H. Austin.  Statics and dynamics of single dna molecules confined in nanochannels.  Phys. Rev. Lett., 94:196101, 2005.
  * [21] T. Odijk.  On the statistics and dynamics of confined or entangled stiff polymers.  Macromolecules, 16:1340, 1983.
  * [22] P Pincus.  Excluded volume effects and stretched polymer chains.  Macromolecules, 9:386–388, 1976.
  * [23] Y. Yang, T Burkhardt, and G. Gompper.  Free energy and extension of a semiflexible polymer in cylindrical confining geometries.  Phys. Rev. E, 76:011804, 2007.
  * [24] M. Dijkstra, D. Frenkel, and H. NW Lekkerkerker.  Confinement free energy of semiflexible polymers.  Physica A, 193:374, 1993.
  * [25] K. Jo, D.M Dalia, T Odijk, J. J. de Pablo, M. D. Graham, R Runnheim, D Forrest, and D. C. Schwartz.  A single-molecule barcoding system using nanoslits for DNA analysis.  Proc. Natl. Acad. Sci., 104:2673, 2007.
  * [26] M. Doi and S. F. Edwards.  The Theory of Polymer Dynamics.  Clarendon Press, Oxford, 1988.
  * [27] G. Morrison and D. Thirumalai.  Semiflexible chains in confine spaces.  Phys. Rev. E, 79:011924, 2009.
  * [28] BY Ha and D Thirumalai.  Semiflexible chains under tension.  J. Chem. Phys., 106:4243, 1997.
  * [29] J. Wang and H. Gao.  Stretching a stiff polymer in a tube.  J. Mater. Sci., 42:8838–8843, 2007.
  * [30] Greg Morrison.  Correlation functions for strongly confined wormlike chains.  https://arxiv.org/abs/2003.12187, 2023.
  * [31] BY Ha and D. Thirumalai.  A mean-field model for semiflexible chains.  J. Chem. Phys., 103:9408, 1995.
  * [32] M Hinczewski and RR Netz.  Anisotropic hydrodynamic mean-field theory for semiflexible polymers under tension.  Macromol, 44:6972, 2011.
  * [33] R. P. Feynman.  Statistical Mechanics.  W. A. Benjamin, Inc., 1982.