# Calorimetry of a Quantum Phase Slip

E. Gümüş  Univ. Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 25 rue des
Martyrs, Grenoble, France  D. Majidi  Univ. Grenoble Alpes, CNRS, Grenoble
INP, Institut Néel, 25 rue des Martyrs, Grenoble, France  D. Nikolić
Fachbereich Physik, Universität Konstanz, D-78457 Konstanz, Germany  P. Raif
Univ. Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 25 rue des Martyrs,
Grenoble, France  Fachbereich Physik, Universität Konstanz, D-78457 Konstanz,
Germany  B. Karimi  QTF Centre of Excellence, Department of Applied Physics,
Aalto University School of Science, P.O. Box 13500, 00076 Aalto, Finland  J.
T. Peltonen  QTF Centre of Excellence, Department of Applied Physics, Aalto
University School of Science, P.O. Box 13500, 00076 Aalto, Finland  E. Scheer
Fachbereich Physik, Universität Konstanz, D-78457 Konstanz, Germany  J. P.
Pekola  QTF Centre of Excellence, Department of Applied Physics, Aalto
University School of Science, P.O. Box 13500, 00076 Aalto, Finland  H.
Courtois  Univ. Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 25 rue des
Martyrs, Grenoble, France  W. Belzig  Fachbereich Physik, Universität
Konstanz, D-78457 Konstanz, Germany  C. B. Winkelmann  Univ. Grenoble Alpes,
CNRS, Grenoble INP, Institut Néel, 25 rue des Martyrs, Grenoble, France

###### Abstract

In a Josephson junction, which is the central element in superconducting
quantum technology, irreversibility arises from abrupt slips of the gauge-
invariant quantum phase difference across the contact. A quantum phase slip
(QPS) is often visualized as the tunneling of a flux quantum in the transverse
direction to the superconducting weak link, which produces dissipation. In
this work, we detect the instantaneous heat release caused by a QPS in a
Josephson junction using time-resolved electron thermometry on a
nanocalorimeter, signaled by an abrupt increase of the local electronic
temperature in the weak link and subsequent relaxation back to equilibrium.
Beyond providing a cornerstone in experimental quantum thermodynamics in form
of observation of heat in an elementary quantum process, this result sets the
ground for experimentally addressing the ubiquity of dissipation, including
that in superconducting quantum sensors and qubits.

The magnetic flux threading a superconducting loop is quantized in units of
the flux quantum $\Phi_{0}=h/2e$. The tunneling of a flux quantum in or out of
such a loop is associated to a change of $2\pi$ in the winding of the phase of
the quantum wave function along the loop. The manipulation of individual flux
quanta is at the core of superconducting circuit logics, both in the classical
and in the quantum information regime Mooij and Harmans (2005); Mooij and
Nazarov (2006); Liebermann and Wilhelm (2016). Rapid single-flux-quantum
(RSFQ) logic can operate up to 100 GHz frequencies and is considered as
promising classical control electronics of qubits Leonard  _et al._ (2019);
Howington  _et al._ (2019); McDermott  _et al._ (2018). In the quantum regime,
the coherent superposition and manipulation of flux states is at the basis of
flux qubits Chiorescu  _et al._ (2003); Yan  _et al._ (2016) and the fluxonium
Manucharyan  _et al._ (2009); Nguyen  _et al._ (2019).

Slips of the quantum phase occur when the gauge-invariant phase difference
across a weak link, that is, a Josephson junction, in the superconducting loop
suddenly relaxes (Fig. 1). Quantum phase slips (QPSs) are ubiquitous in
superconducting electronics and can be seen as the dual process to Cooper pair
tunneling. Furthermore, coherent QPSs have been proposed as a building block
for phase-slip qubit devices Mooij and Harmans (2005); Mooij and Nazarov
(2006); Astafiev  _et al._ (2012). On the other hand, their proliferation is
responsible for the destruction of superconductivity in one dimension Golubev
and Zaikin (2001) and can lead to thermal avalanches in current-biased
Josephson junctions Sahu  _et al._ (2009). In essence, a QPS can be considered
the quantum of dissipation in superconducting electronics.

In this work, we investigate the thermal signature of incoherent individual
quantum phase slips in a superconducting Josephson junction. Applying a MHz-
bandwidth electron thermometry technique to a Josephson junction, we measure
the heat generated by a single QPS as well as the subsequent thermal
relaxation. The data are in good agreement with a theoretical model that we
developed for describing the superconducting properties of the device. Our
work therefore demonstrates the possibility to quantitatively account for
dissipative effects in quantum nanoelectronics, with evident applications to
quantum computing.

Figure 1: Quantum phase slip in a Josephson junction. (a) Real-space sketch of
the QPS mechanism: at the instability point of the $\Phi_{x}(\varphi)$
relation, the phase drop $\varphi$ and the screening current $I_{s}$ relax
abruptly to smaller values, as a quantum of flux tunnels perpendicular to the
Josephson junction (dark grey), releasing heat. (b) Phase drop $\varphi$
across the SNS junction versus applied flux to the SQUIPT, following Eq. (1)
with $\beta=10$. The dashed part of the curve cannot be accessed. In a quantum
phase slip (blue arrows), $\varphi$ changes by slightly less than $2\pi$. (c)
Potential energy of the SQUIPT as a function of $\varphi$. A local energy
minimum can become unstable as the externally applied flux is changed. By
macroscopic quantum tunneling of the phase, a lower energy valley is reached,
releasing an energy $\Delta U$.

The experimental core element is a superconducting quantum interference device
(SQUID), that is, a superconducting loop containing one (as is the case here)
or more Josephson junctions, and to which a magnetic flux $\Phi_{x}$ is
applied. The difference between the applied and the physical flux is absorbed
by screening supercurrents, leading to a gradient in the phase of the quantum
wave function along the loop, and, to a large extent, across the Josephson
weak link. The relation between $\Phi_{x}$ and the phase drop $\varphi$ across
the Josephson junction can be written as

$2\pi(\Phi_{x}/\Phi_{0})=\varphi+\beta\sin(\varphi),$ (1)

with the screening parameter $\beta=2\pi LI_{c}/\Phi_{0}$, where $L$ is the
loop inductance and $I_{c}$ the Josephson junction’s critical current Tinkham
(1996). Irreversibility arises when the SQUID’s magnetic screening parameter
$\beta$ exceeds 1 and Eq. (1) is no longer single-valued, as illustrated in
Fig. 1b,c. In this situation, the penetration of an additional flux quantum
into the SQUID loop does not occur smoothly and reversibly, but via a sudden
tunneling of the phase, that is, a quantum phase slip.

Figure 2: Hysteretic RF-SQUIPT. (a) Device schematics, including a false-color
scanning electron micrograph (scale bar 5 $\mu$m) of the SQUIPT loop (yellow
area) and the rapid flux line (pink). The zoom highlights the SNS junction (Al
= blue, Cu = grey), connected laterally by the tunnel contact (cyan) (scale
bar 200 nm), connected to the resonator. (b) Bias spectroscopy $P_{\rm
out}(V_{b})$ at cryostat temperatures 50 mK (blue) and 400 mK (red). (c)
Calibration of $P_{\rm out}(V_{b}=0)$ versus cryostat temperature at
equilibrium, under two different phase drops. The grey shaded region thus
covers all possible values of $P_{\rm out}(T,\varphi)$. (d) Resonator response
at $V_{b}=0$ as a function of increasing applied magnetic flux, at three
cryostat temperatures (50 mK, blue; 200 mK, black; 400 mK, red). The dark red
line exemplifies the response to a downward sweep of the flux. The dashed line
is a calculation (see text). (e) Temperature dependence of the screening
parameter $\beta$, extracted from (d), and theoretical fit (red line, see
text).

The above concepts are the basis of extensive applications in superconducting
electronics, such as SQUID magnetometry and superconducting flux qubits van
Der Wal  _et al._ (2000); Friedman  _et al._ (2000); Chiorescu  _et al._
(2003). Notably, the Josephson weak links used in such SQUIDs can be provided
by a variety of junction types, including tunnel junctions, micro-bridges, and
proximity weak links Likharev (1979). In the latter, a short, non-
intrinsically-superconducting element, such as a normal metallic wire allows
for superconducting correlations to propagate between both superconducting
reservoirs. Due to the wire’s usually low normal state resistance, following
the resistively and capacitively shunted Josephson junction model, the quantum
phase dynamics in such superconductor-normal metal-superconductor (SNS)
Josephson junctions are inherently overdamped Tinkham (1996). This ensures
that upon a QPS, the quantum phase only evolves towards the nearest
neighboring potential valley (Fig. 1c). This releases an energy

$\Delta U=\frac{\Phi_{0}}{2\pi L}\int\Phi_{x}\,d\varphi,$ (2)

which depends on the magnitude of the quantum phase slip. In the large $\beta$
limit, the phase jumps by about $2\pi$, and $\Delta U\approx I_{c}\Phi_{0}$.

Our device consists of a SQUID with a single Josephson junction, provided by
an SNS weak link of length $520$ nm (Fig. 2a). Here, the superconducting
circuit parts are made of aluminum, while the metal N is made of copper. While
the loop is grounded, the center of N is further contacted by another
superconducting finger through a tunnel junction (Fig. 2a), with normal state
resistance $R_{T}\approx 7$ k$\Omega$. This SQUID variant was named SQUIPT
Giazotto  _et al._ (2010), where PT stands for proximity transistor. Our
SQUIPT was designed to be in the hysteretic regime. For this, we took into
account the geometric and kinetic inductance contributions, leading here to
$L\approx 630\pm 50$ pH. The SNS junction critical current $I_{c}$ cannot be
determined independently in this device, but is expected from similar SNS
devices in a current-biased geometry to be of a few $\mu$A Courtois  _et al._
(2008); Angers  _et al._ (2008); Dutta  _et al._ (2020). From this, values of
the screening parameter $\beta\sim 10$ can be anticipated, in good agreement
with experiments, as discussed below.

Applying a dc voltage $V_{b}$ to the tunnel junction, we can perform tunnel
spectroscopy by measuring its differential conductance $G(V_{b})$. Here, $G$
is not read out, as usual, by a low-frequency transport measurement, but by a
radio-frequency (RF) technique, using a superconducting $LC$ resonator with
resonant frequency 575 MHz Schmidt  _et al._ (2003); Gasparinetti  _et al._
(2015); Karimi and Pekola (2018). By embedding the tunnel junction in parallel
to the resonator, and for a fixed incident RF power $P_{\rm in}$, changes in
$G$ can be read out by their effect on the transmitted power at resonance
$P_{\rm out}$, which we record after cryogenic amplification (details in Supp.
Info. file). This has the paramount advantage of allowing for extremely rapid
measurements, limited by the resonator bandwidth, of about 10 MHz here.

Figure 2b shows measurements of $P_{\rm out}(V_{b})$ at two temperatures.
Several characteristic spectroscopic features stand out, in particular (i) a
spectroscopic gap of total width 480 $\mu$V, (ii) subgap resonances near $\pm
190\,\mu$V visible only at 400 mK, and (iii) three low-energy resonances at 0
and $\pm 15\,\mu$V, respectively. Keeping in mind that the tunnel junction
connects an intrinsic superconductor with gap $\Delta$ and a proximized metal
with a (smaller) induced gap $E_{g}$, we can evaluate the total spectroscopic
gap as $2(\Delta+E_{g})/e$. At intermediate temperatures, thermally activated
conductance resonances occur at $\pm|\Delta-E_{g}|/e$, which we identify as
feature (ii). From the two above relations we find a gap $\Delta=210\pm
5\,\mu$eV in the superconducting probe electrode at 400 mK (225 $\mu$eV at 50
mK), a typical value in nanostructured aluminum. We further extract
$E_{g}=29\,\mu$eV, in good agreement with theoretical estimates (details in
the Supp. Info. file). Eventually, the conductance resonance at $V_{b}=0$
(feature (iii)) is a signature of the Josephson coupling across the tunnel
junction, which was purposely designed to have an intermediate transparency
Karimi and Pekola (2018). The satellite peaks at $\pm 15\,\mu$V are probably
the signature of inelastic Cooper pair tunneling Karimi and Pekola (2018) and
are not of central relevance to this work.

As visible in Fig. 2b,c and already discussed in detail in Karimi and Pekola
(2018); Karimi  _et al._ (2020a, b), the zero-bias conductance of the tunnel
junction, and therefore $P_{\rm out}(V_{b}=0)$, is a sensitive probe of the
electron temperature $T$ in N. Accordingly, we set $V_{b}=0$ in the remainder
of this work and use $P_{\rm out}$ for both static and dynamic electron
thermometry, after initial calibration under equilibrium conditions (Fig. 2c).
However, and in contrast to previous work Karimi and Pekola (2018); Karimi
_et al._ (2020a, b) in which N was not subject to a phase drop, in the present
device the tunnel junction conductance and thus $P_{\rm out}$ are clearly also
a function of the phase $\varphi$ across the SNS junction. This is seen in
Fig. 2d, where sweeping the applied flux translates into a phase variation via
Eq. (1). As $\varphi\to\pi/2$, the decrease of $E_{g}$ Zhou  _et al._ (1998);
Le Sueur  _et al._ (2008) entails a rapidly shrinking Josephson energy $E_{J}$
of the tunnel junction, and thus a decrease of $G(0)$, that is, an increase of
$P_{\rm out}$. When the switching point in the $\Phi_{x}(\varphi)$ relation is
reached, $\varphi$ suddenly relaxes to a smaller value (modulo $2\pi$),
restoring $E_{J}$ and thus leading to an abrupt drop in $P_{\rm out}$. As
expected, the same pattern is repeated with period $\Phi_{0}$ in the applied
flux and mirror symmetric under inversion of the sweep direction Ligato  _et
al._ (2021). At higher temperatures, $\beta\propto I_{c}$ decreases and the
modulation amplitude of $P_{\rm out}(\Phi_{x})$ shrinks, while the
$\Phi_{0}$-periodicity of the signal is preserved.

Figure 3: Heat relaxation dynamics after a QPS. (a) Time-resolved electron
temperature in the absorber, at different starting temperatures set by the
cryostat bath, following a $70$-ns flux pulse at $t=0$. The data sampling rate
is 2 MHz, the data shown are the result of averaging over $10^{5}$ pulses, at
a repetition rate of 2 kHz. (b) Return to equilibrium $\Delta T(t)$ at 100 mK,
following a flux pulse. Same data as the black curve in (a), but in a semi-
log-scale representation and over a wider time window. The red line is a
calculation based on the model discussed in the text. The dashed line is an
exponential fit (with time constant $\tau\approx 140\,\mu$s) to the long-term
relaxation, evidencing the presence of a second slowly relaxing bath. (c)
Magnitude of the initial temperature rise $\Delta T_{0}$ at $t=0$, determined
by the fit as shown in (b) (bullets). The solid line is a calculation (see
text).

For a quantitative understanding of the RF-SQUIPT, we use the quasi-classical
Usadel equations Larkin and Ovchinnikov (1968); Usadel (1970); Zhou  _et al._
(1998); Belzig  _et al._ (1999), with a single consistent set of microscopic
parameters, described in detail in the Supp. Info. file. The density of states
in N is known to display a minigap, which depends approximately on $\varphi$
like $E_{g}(\varphi)=E_{g}(0)\,|\cos(\varphi/2)|$ Zhou  _et al._ (1998). The
tunnel junction connecting the condensate in N to the superconducting probe
electrode has a Josephson energy $E_{J}(\varphi,T)$ and thus a zero-bias
conductance $G(\varphi,T)$, which can be drawn back analytically to
$E_{g}(\varphi)$ and the critical current $I_{c}$ of the SNS junction. In
combination with Eq. (1) and the relation between $G$ and $P_{\rm out}$, the
calculation provides an accurate description of the applied-flux dependence of
the RF signal (Fig. 2d), the only adjustable parameter being the magnitude of
$\beta$.

The temperature dependence of $\beta$ extracted from the data in Fig. 2d is
plotted in Fig. 2e, following the trend expected for $I_{c}(T)$ in an SNS
junction. The solid line shows the calculation from the same model as above,
yielding a $5.9\,\mu$A zero-temperature critical current in the SNS junction.
The parameters entering the calculation, and in particular the minigap, are
determined independently using the tunnel spectra (Fig. 2b). We attribute the
slight underestimation of $\beta$ by theory above 300 mK to the temperature
dependence of the kinetic inductance, which was not accounted for.

The data discussed so far, and summarized in Fig. 2, provide a consistent
physical understanding of the RF-SQUIPT under time-averaged, and thus
isothermal conditions. We now move to the time-resolved response $P_{\rm
out}(t)$, which displays the calorimetric signature of the heat deposited by
individual QPS events. In addition to the static flux bias, we apply a time-
dependent (square-wave or pulsed) flux, by passing a current $I_{\rm RFL}(t)$
through the superconducting rapid flux line (RFL), visible on the left side of
the SQUIPT in Fig. 2a. The instantaneous flux bias is $\Phi_{x}(t)=M\,I_{\rm
RFL}(t)$, with a mutual inductance $M=12.1$ pH further discussed in the Supp.
Info. file. In order to increase the signal-to-noise ratio, we average the
resonator response over a large number of periodically generated identical
$70$ ns flux pulses.

As long as the amplitude of the flux pulses does not exceed the threshold
leading to a QPS, $P_{\rm out}$ follows changes in $\Phi_{x}(t)$
instantaneously (not shown, see Supp. Info. file for details). However, as
soon as the threshold to instability of the flux state is overcome, $P_{\rm
out}(t)$ displays in addition a slower relaxation to its novel equilibrium,
which is indicative of thermalization. As evidenced from Fig. 3a, the
relaxation dynamics after a flux pulse are strongly temperature dependent, as
expected for instance from a dominantly electron-phonon (e-ph) coupling-driven
thermalization after an initial heating event Gasparinetti  _et al._ (2015).
Above 300 mK, the relaxation times become too short ($<1\,\mu$s) to be
measured.

The thermal dynamics can be described by a heat balance equation, basing on
the standard assumption that variations of the absorber’s internal energy $U$
are evacuated to a heat bath via electron-phonon coupling. In metallic
nanostructures at low temperatures, this power is usually written
${\dot{Q}}_{e-ph}=\Sigma{\mathcal{V}}(T^{5}-T_{0}^{5})$, where $\Sigma=2\times
10^{9}$ W K-5 m-3 is the electron-phonon coupling constant in copper Giazotto
_et al._ (2006); Viisanen and Pekola (2018), and ${\mathcal{V}}\approx 8\times
10^{-21}$ m-3 the geometrically estimated absorber volume. Given the rather
small $E_{g}$ in the SNS junction, our calculations indicate that the
proximity effect should only lead to negligibly small departures from the
normal-state electron-phonon coupling and heat capacity $C$, at the
experimental temperatures. Therefore we write $U=\gamma T^{2}/2$, using a
reported $\gamma=71$ J m-3 K-2 in nanostructured Cu Viisanen and Pekola
(2018). As seen in Fig. 3a, the temperature increase after the initial heat
pulse can be large compared to the starting temperature. The dynamical heat
balance differential equation thus cannot be linearized in $\Delta T=T-T_{b}$,
and must be solved numerically. The result at $T_{b}=100$ mK is shown in Fig.
3b (red line), where the only free parameter is the initial temperature
increase $\Delta T_{0}$. The calculation follows closely the data during the
first period of the thermal relaxation. Interestingly, the initial fast decay
is rapidly taken over by a much slower process, which was already reported for
nanoscale Cu absorbers Viisanen and Pekola (2018); Wang  _et al._ (2019). This
might be due to another heat reservoir, for instance in surface states of the
Cu absorber.

The initial temperature increase $\Delta T_{0}$ after a flux pulse is plotted
as a function of the starting temperature in Fig. 3c. Naturally, the
temperature rise is highest when starting from the lowest base temperature,
because at higher temperature both $C$ increases and $\Delta U\sim
I_{c}\Phi_{0}$ decreases. For a quantitative modeling, one must bear in mind
that a flux pulse induces necessarily two QPSs: one while ramping up the flux
bias, and the second when returning to the initial state. We account for the
fact that the return QPS occurs at a higher electronic temperature and thus
produces less dissipation. The sudden temperature rise after the first QPS
could also lead the sinusoidal potential energy landscape modulation (see Fig.
1c) to collapse to the point that the reached valley also becomes unstable and
the phase eventually jumps by $4\pi$, or more. This was actually observed in
another device, which we describe in the Supp. Info. file, but does not apply
in the sample described here. Finally, the quasi-classical Usadel formalism
with the same set of microscopic parameters as previously, in combination with
Eqs. (1,2) and the known values of $\gamma$ and ${\mathcal{V}}$, describes
very accurately (red line in Fig. 3c) the initial temperature rise $\Delta
T_{0}$.

Our findings of microsecond-scale thermal relaxation following a QPS tracked
in real time highlight both the paramount effect dissipation can produce in
quantum circuits, and the potential of large-bandwidth electron thermometry
for quantum thermodynamics in nanoelectronic circuits. The detection of a
quantum of dissipation in a Josephson junction opens promising perspectives
for future experiments. These would aim, for instance, at detecting the minute
dissipation arising from a projective qubit measurement with emission of a
single microwave photon Pekola and Karimi (2022), or the detection of the
elusive axion with even smaller energy Braine  _et al._ (2020). On the other
hand, the proper measurement and control of dissipation in QPS-based circuits
will allow for overcoming self-heating-limited device performances.

We acknowledge help from A. Théry and T. Crozes. The samples were fabricated
at the Nanofab platform at Institut Néel. This work received support from the
European Union under the Marie Sklodowska-Curie Grant Agreement No. 766025
(QuESTech), from the Agence Nationale de la Recherche under the program
“Investissements d’avenir” (ANR-15-IDEX-02), from the Laboratoire d’excellence
LANEF (ANR-10-LABX-51-01), and from Deutsche Forschungsgemeinschaft (DFG;
German Research Foundation) via SFB 1432 (Project No. 425217212).

## References

  * Mooij and Harmans (2005) J. E. Mooij and C. J. P. M. Harmans, “Phase-slip flux qubits,” New J. Phys. 7, 219 (2005).
  * Mooij and Nazarov (2006) J. E. Mooij and Yu. V. Nazarov, “Superconducting nanowires as quantum phase-slip junctions,” Nat. Phys. 2, 169–172 (2006).
  * Liebermann and Wilhelm (2016) P. J. Liebermann and F. K. Wilhelm, “Optimal qubit control using single-flux quantum pulses,” Phys. Rev. Appl. 6, 024022 (2016).
  * Leonard  _et al._ (2019) E. Leonard, M. A. Beck, J. Nelson, B. G. Christensen, T. Thorbeck, C. Howington, A. Opremcak, I. V. Pechenezhskiy, K. Dodge, N. P. Dupuis, M. D. Hutchings, J. Ku, F. Schlenker, J. Suttle, C. Wilen, S. Zhu, M. G. Vavilov, B. L. T. Plourde, and R. McDermott, “Digital Coherent Control of a Superconducting Qubit,” Phys. Rev. Appl. 11, 014009 (2019).
  * Howington  _et al._ (2019) C. Howington, A. Opremcak, R. McDermott, A. Kirichenko, O. A. Mukhanov, and B. L. T. Plourde, “Interfacing Superconducting Qubits With Cryogenic Logic: Readout,” IEEE Transactions on Applied Superconductivity 29, 1–5 (2019).
  * McDermott  _et al._ (2018) R. McDermott, M. G. Vavilov, B. L. T. Plourde, F. K. Wilhelm, P. J. Liebermann, O. A. Mukhanov, and T. A. Ohki, “Quantum–classical interface based on single flux quantum digital logic,” Quantum Science and Technology 3, 024004 (2018).
  * Chiorescu  _et al._ (2003) I. Chiorescu, Y. Nakamura, C. J. P. Harmans, and J. E. Mooij, “Coherent quantum dynamics of a superconducting flux qubit,” Science 299, 1869–1871 (2003).
  * Yan  _et al._ (2016) F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, _et al._ , “The flux qubit revisited to enhance coherence and reproducibility,” Nat. Commun. 7, 1–9 (2016).
  * Manucharyan  _et al._ (2009) V. E. Manucharyan, J. Koch, L. I. Glazman, and M. H. Devoret, “Fluxonium: Single Cooper-pair circuit free of charge offsets,” Science 326, 113 (2009).
  * Nguyen  _et al._ (2019) L. B. Nguyen, Y.-H. Lin, A. Somoroff, R. Mencia, N. Grabon, and V. E Manucharyan, “High-coherence fluxonium qubit,” Phys. Rev. X 9, 041041 (2019).
  * Astafiev  _et al._ (2012) O. V. Astafiev, L. B. Ioffe, S. Kafanov, Y. A. Pashkin, K. Y. Arutyunov, D. Shahar, O. Cohen, and J. S. Tsai, “Coherent quantum phase slip,” Nature 484, 355–358 (2012).
  * Golubev and Zaikin (2001) D. S. Golubev and A. D. Zaikin, “Quantum tunneling of the order parameter in superconducting nanowires,” Phys. Rev. B 64, 014504 (2001).
  * Sahu  _et al._ (2009) M. Sahu, M. Bae, A. Rogachev, D. Pekker, T.-C. Wei, N. Shah, P. M. Goldbart, and A. Bezryadin, “Individual topological tunnelling events of a quantum field probed through their macroscopic consequences,” Nat. Phys. 5, 503–508 (2009).
  * Tinkham (1996) M. Tinkham,  _Introduction to Superconductivity_ (Dover, New York, 2nd edition, 1996).
  * van Der Wal  _et al._ (2000) C. H. van Der Wal, A. C. J. Ter Haar, F. K. Wilhelm, R. N. Schouten, C. J. P. M. Harmans, T. P. Orlando, S. Lloyd, and J. E. Mooij, “Quantum superposition of macroscopic persistent-current states,” Science 290, 773–777 (2000).
  * Friedman  _et al._ (2000) J. R. Friedman, V. Patel, W. Chen, S. K. Tolpygo, and J. E. Lukens, “Quantum superposition of distinct macroscopic states,” Nature 406, 43–46 (2000).
  * Likharev (1979) K. K. Likharev, “Superconducting weak links,” Rev. Mod. Phys. 51, 101 (1979).
  * Giazotto  _et al._ (2010) F. Giazotto, J. T. Peltonen, M. Meschke, and J. P Pekola, “Superconducting quantum interference proximity transistor,” Nat. Phys. 6, 254–259 (2010).
  * Courtois  _et al._ (2008) H. Courtois, M. Meschke, J. T. Peltonen, and J. P Pekola, “Origin of hysteresis in a proximity Josephson junction,” Phys. Rev. Lett. 101, 067002 (2008).
  * Angers  _et al._ (2008) L. Angers, F. Chiodi, G. Montambaux, M. Ferrier, S. Guéron, H. Bouchiat, and J. C. Cuevas, “Proximity dc squids in the long-junction limit,” Phys. Rev. B 77, 165408 (2008).
  * Dutta  _et al._ (2020) B. Dutta, D. Majidi, N. W. Talarico, N. Lo Gullo, H. Courtois, and C. B. Winkelmann, “Single-quantum-dot heat valve,” Phys. Rev. Lett. 125, 237701 (2020).
  * Schmidt  _et al._ (2003) D. R. Schmidt, C. S. Yung, and A. N. Cleland, “Nanoscale radio-frequency thermometry,” Appl. Phys. Lett. 83, 1002–1004 (2003).
  * Gasparinetti  _et al._ (2015) S. Gasparinetti, K. L. Viisanen, O.-P. Saira, T. Faivre, M. Arzeo, M. Meschke, and J. P. Pekola, “Fast electron thermometry for ultrasensitive calorimetric detection,” Phys. Rev. Appl. 3, 014007 (2015).
  * Karimi and Pekola (2018) B. Karimi and J. P Pekola, “Noninvasive thermometer based on the zero-bias anomaly of a superconducting junction for ultrasensitive calorimetry,” Phys. Rev. Appl. 10, 054048 (2018).
  * Karimi  _et al._ (2020a) B. Karimi, D. Nikolić, T. Tuukkanen, J. T. Peltonen, W. Belzig, and J. P Pekola, “Optimized proximity thermometer for ultrasensitive detection,” Phys. Rev. Appl. 13, 054001 (2020a).
  * Karimi  _et al._ (2020b) B. Karimi, F. Brange, P. Samuelsson, and J. P Pekola, “Reaching the ultimate energy resolution of a quantum detector,” Nat. Commun. 11, 1–6 (2020b).
  * Zhou  _et al._ (1998) F. Zhou, P. Charlat, B. Spivak, and B. Pannetier, “Density of states in superconductor-normal metal-superconductor junctions,” J. Low Temp. Phys. 110, 841–850 (1998).
  * Le Sueur  _et al._ (2008) H. Le Sueur, P. Joyez, H. Pothier, C. Urbina, and D. Estève, “Phase controlled superconducting proximity effect probed by tunneling spectroscopy,” Phys. Rev. Lett. 100, 197002 (2008).
  * Ligato  _et al._ (2021) N. Ligato, E. Strambini, F. Paolucci, and F. Giazotto, “Preliminary demonstration of a persistent Josephson phase-slip memory cell with topological protection,” Nat. Commun. 12, 1–8 (2021).
  * Larkin and Ovchinnikov (1968) A. L. Larkin and Yu. N. Ovchinnikov, “Quasiclassical method in the theory of superconductivity,” Zh. Eksp. Teor. Fiz. 55, 2262 (1968).
  * Usadel (1970) K. D. Usadel, “Generalized diffusion equation for superconducting alloys,” Phys. Rev. Lett. 25, 507 (1970).
  * Belzig  _et al._ (1999) W. Belzig, F. K. Wilhelm, C. Bruder, G. Schön, and A. D. Zaikin, “Quasiclassical Green’s function approach to mesoscopic superconductivity,” Superlattices Microstruct. 25, 1251 (1999).
  * Giazotto  _et al._ (2006) F. Giazotto, T. T. Heikkilä, A. Luukanen, A. M. Savin, and J. P. Pekola, “Opportunities for mesoscopics in thermometry and refrigeration: Physics and applications,” Rev. Mod. Phys. 78, 217 (2006).
  * Viisanen and Pekola (2018) K. L. Viisanen and J. P. Pekola, “Anomalous electronic heat capacity of copper nanowires at sub-kelvin temperatures,” Phys. Rev. B 97, 115422 (2018).
  * Wang  _et al._ (2019) L. B. Wang, D. S. Golubev, Y. M. Galperin, and J. P. Pekola, “Dynamic thermal relaxation in metallic films at sub-kelvin temperatures,” arXiv preprint arXiv:1910.09448 (2019).
  * Pekola and Karimi (2022) J. P. Pekola and B. Karimi, “Ultrasensitive calorimetric detection of single photons from qubit decay,” Phys. Rev. X 12, 011026 (2022).
  * Braine  _et al._ (2020) T. Braine, R. Cervantes, N. Crisosto, N. Du, S. Kimes, L. J. Rosenberg, G. Rybka, J. Yang, D. Bowring, A. S. Chou, _et al._ , “Extended search for the invisible axion with the axion dark matter experiment,” Phys. Rev. Lett. 124, 101303 (2020).


# Modified (hyper-)viscosity for coarse-resolution ocean models

Louis Thiry, Long Li and Etienne Mémin
INRIA/IRMAR Rennes, France

###### Abstract

We present a simple parameterization for coarse-resolution ocean models. To
replace computationally expensive high-resolution ocean models, we develop a
computationally cheap parameterization for coarse-resolution models based
solely on the modification of the viscosity term in advection equations. It is
meant to reproduce the mean quantities like pressure, velocity, or vorticity
computed from a high-resolution reference solution or using observations. We
test this new parameterization on a double-gyre quasi-geostrophic model in the
eddy-permitting regime. Our results show that the proposed scheme improves
significantly the energy statistics and the intrinsic variability on the
coarse mesh. This method shall serve as a deterministic basis model for
coarse-resolution stochastic parameterizations in future works.

##  1 Introduction

Ocean general circulation models used at climatic scales are limited for
evident computational reasons to too coarse horizontal resolutions to solve
correctly ocean mesoscale and sub-mesoscale eddies, even with large
computational infrastructure. The horizontal resolution of the most recent
climatic ocean models is of the order of the Rossby radius of deformation.
These models are hence in the so-called eddy-permitting regime and they can
solve partially the mesoscale (i.e. 10-100 km) eddy field. These models
however suffer from strong limitations. In particular, they are unable to
reproduce accurately large-scale structures such as the eastward turbulent jet
in an idealized double-gyre configuration.

Recent parameterizations have shown significant improvements in coarse-
resolution models compared to high-resolution reference solutions [3].
However, it remains an important topic of research, as the actual generation
of parametrizations is not completely able to resolve the effects of the
unresolved scales on the large-scale flow structures.

A wide range of subgrid parametrizations relies on eddy viscosity such as
Laplacian and biharmonic schemes [17, 11, 5, 4]. It has been shown in [10]
that including only these (hyper)viscosity in coarse-resolution models often
causes too much dissipation and results in an artificial energy sink at large
scales. In general, even eddy-permitting models are not energetic enough and
as a result, the long-time average of any coarse model’s variable of interest
departs completely from the long-time average of high-resolution models
subsampled at the same scale. This becomes the main motivation of the present
work. In particular, we would like to answer the following question: how can
we reduce the excessive resolved kinetic energy loss due to the viscosity
while simultaneously ensuring numerical stability?

We propose a simple affine parameterization of (hyper)viscosity. The
(bi)laplacian operator $\Delta^{p}f$ is replaced by
$\Delta^{p}\left(f-f^{\prime}\right)$, where $f^{\prime}$ is a field of same
dimension as $f$ that does not depend upon time. We interpret this method as a
mathematical regularization technique to guide the solutions towards prior
information. We frame $f^{\prime}$ as the solution of an optimal control
problem to reproduce statistics computed from a reference solution or
observations. We present a method to solve this optimal control problem.

We test the proposed method with an idealized double-gyre configuration. For
that purpose, we release with this article a fast, concise, and CPU-GPU
portable Pytorch implementation of a multi-layer quasi-geostrophic model on a
rectangular domain. We implement and test our optimization procedure within
this setting.

This article is organized as follows: we present in section 2 the double gyre
quasi-geostrophic model we use and detail its implementation, we present in
section 3 our modified viscosity parameterization and we show and discuss
numerical results in section 4.

##  2 Double gyre quasi-geostrophic model

###  2.1 Governing equations

We use the same multi-layer quasi-geostrophic model in a non-periodic
rectangular domain as in [7]. Here, we only give a brief review of this
system. The quasi-geostrophic pressure and potential vorticity (PV) are
stacked in three isopycnal layers. We adopt vector forms to denote the layered
pressure and potential vorticity (PV):

$\mathbf{p}=\begin{bmatrix}p_{1}\\\ p_{2}\\\
p_{3}\end{bmatrix},\quad\mathbf{q}=\begin{bmatrix}q_{1}\\\ q_{2}\\\
q_{3}\end{bmatrix}\ .$

The forced and damped quasi-geostrophic (QG) equations can be then written as

$\displaystyle\partial_{t}\mathbf{q}=\frac{1}{f_{0}}J(\mathbf{q},\mathbf{p})+f_{0}B\mathbf{e}+\frac{1}{f_{0}}\big{(}a_{2}\Delta-
a_{4}\Delta^{2}\big{)}(\Delta\mathbf{p}),$ (1) $\displaystyle\big{(}\Delta-
f_{0}^{2}A\big{)}\mathbf{p}=f_{0}\mathbf{q}-f_{0}\beta(y-y_{0}),$ (2)

where $\Delta=\partial_{xx}^{2}+\partial_{yy}^{2}$ denotes the horizontal
Laplacian, $\Delta^{2}$ the bi-laplacian operator,
$J(a,b)=\partial_{x}a\partial_{y}b-\partial_{x}b\partial_{y}a$ stands for the
Jacobi operator, $f_{0}+\beta(y-y_{0})$ is the Coriolis parameter under beta-
plane approximation with the meridional axis center $y_{0}$, $a_{2}$ and
$a_{4}$ are the Laplacian and biharmonic viscosity coefficients. Parameters of
the configuration are listed in the Tables 1 and 2 in Appendix. Besides, the
second term on the right-hand side of equation ((1)) represents the external
forcing applied on different layers. In this work, we only consider an
idealized case in which the ocean basin is driven by a stationary and
symmetric wind stress $\vec{\tau}=(\tau^{x},\tau^{y})$ on the surface and by a
linear Ekman stress at the bottom. In that case, the forcing term can be
specified by

$B=\begin{bmatrix}\frac{1}{H_{1}}&\frac{-1}{H_{1}}&0&0\\\
0&\frac{1}{H_{2}}&\frac{-1}{H_{2}}&0\\\
0&0&\frac{1}{H_{3}}&\frac{-1}{H_{3}}\\\ \end{bmatrix},\
\quad\mathbf{e}=\begin{bmatrix}\partial_{x}\tau^{y}-\partial_{y}\tau^{x}\\\
0\\\ 0\\\ \frac{\delta_{\rm ek}}{2|f_{0}|}\Delta p_{3}\end{bmatrix},\
\quad\vec{\tau}=\tau_{0}\begin{bmatrix}-\cos(2\pi y/L_{y})\\\ 0\end{bmatrix},$

where $\tau_{0}$ is the magnitude of surface wind, $H_{k}$ is the background
thickness of layer $k$, and $\delta_{\rm ek}$ is the bottom Ekman layer
thickness. The vertical stratification level of such a model is described by
the term $-f_{0}^{2}A\mathbf{p}$ in Equation (2) with

$A=\begin{bmatrix}\frac{1}{H_{1}g_{1.5}^{\prime}}&\frac{-1}{H_{1}g_{1.5}^{\prime}}&0\\\
\frac{-1}{H_{2}g_{1.5}^{\prime}}&\frac{1}{H_{2}}\left(\frac{1}{g_{1.5}^{\prime}}+\frac{1}{g_{2.5}^{\prime}}\right)&\frac{-1}{H_{2}g_{2.5}^{\prime}}\\\
0&\frac{-1}{H_{3}g_{2.5}^{\prime}}&\frac{1}{H_{3}g_{2.5}^{\prime}}\end{bmatrix},$

where $g_{k+0.5}$ is the reduced gravity defined across the interface between
layers $k$ and $k+1$. A multi-layered generalization of this model can be
found in [6]. Note also that such a multi-layered model can be considered as a
vertical discretized approximation of the continuously stratified QG system
[18] with $\partial_{z}(f_{0}\partial_{z}\mathbf{p}/\vec{N}^{2})\approx-
f_{0}A\mathbf{p}$ approximated by finite differences, and in which $\vec{N}$
denotes the buoyancy (or Brunt-Vaisala) frequency.

###  2.2 Pytorch implementation

To facilitate numerical developments and benefit from built-in automatic
differentiation, we develop a Pytorch [13] implementation of the above-
described multilayer QG model 111Available at
https://github.com/louity/qgm_pytorch. For this purpose, we follow rigorously
the strategy of [8]:

  1. 1.

we use a regular numerical grid with finite differences

  2. 2.

We solve the PV advection equation (1) on the whole domain except the
boundaries. We use a standard 5-point finite difference scheme for the
(bi-)laplacian and the energy-enstrophy conservative Arakawa-Lamb scheme for
the Jacobian [1].

  3. 3.

We apply a vertical change of coordinate to equation (2) which becomes a set
of three inhomogeneous Helmholtz equations. We solve these equations with the
spectral Discrete Sine Transform (DST) method, and we add corresponding
homogeneous Helmholtz equation solutions to ensure mass conservation.

  4. 4.

We update the boundary values of the potential vorticity $\mathbf{q}$ using
equation (2).

Detailed equations and numerical routine design choices can be found in [8].
We use a Heun-Runge-Kutta 2 time-stepping instead of the Leap-Frog time scheme
used by [8].

For sake of numerical efficiency, we follow the recommendation of [15]: we
compile computationally demanding routines and simplify finite difference
calculations by reducing as much as possible the number of multiplications. We
end up with a very concise code (less than 300 lines) that only depends upon
Numpy and Pytorch libraries. This implementation will be open-sourced at the
time of the publication.

###  2.3 Eddy-resolving and eddy-permitting regimes

We consider two spatial settings for our simulations :

  1. 1.

The eddy-resolving regime, our high-resolution reference with a 5 km
resolution.

  2. 2.

The eddy-permitting regime, our low-resolution setting with a 40 km
resolution.

Parameters for these two different regimes are written in table 2 in Appendix.

[16] studied the resulting flows’ differences between these two regimes. The
high-resolution eddy-resolving model shows a well-pronounced eastward jet
fuelled by mesoscale eddies circulating while the low-resolution eddy-
permitting model does not induce a proper eastward jet as shown on Fig. 1.
Temporal statistics significantly differ between high- and low-resolution
simulations.

Figure 1: (Top) high-resolution and (bottom) low-resolution top-layer
snapshots after 400 years of integration starting from zero velocity.
Velocities are in m s${}^{-}1$ and PV in s-1.

##  3 Proposed modified viscosity

###  3.1 Motivation

In both resolutions, we use biharmonic viscosity as in [17, 11, 5, 4]
essentially because it is less dissipating at large scales than a Laplacian.
Compared to the usual Laplacian viscosity, it preserves large-scale
structures. However, hyperviscosity remains much too dissipative in the “eddy-
permitting” regime [10]. This too strong dissipation kills the eastward jet
that is present in the high-resolution and that we expect to see in such a
double-gyre quasi-geostrophic model. Fig. 2 shows a sequence of snapshots of
the low-resolution models where we input a downsampled snapshot of the high-
resolution (see Appendix for details on downsampling). After as few as three
years, the eastward jet has almost disappeared, showing that the model is too
dissipating. Lowering the hyper-viscosity coefficient by a factor of 10 does
not solve this problem, and creates spurious gradients in the potential
vorticity as shown in Fig. 2. These numerical artifacts are due to a bad
representation of the direct enstrophy cascade, causing a piling up of the
small-scale vorticity gradients at the cut-off frequency together with
aliasing effects.

Figure 2: (left) Initial condition: high-resolution snapshot on the low-
resolution grid.(center and right) Zonal velocity and potential vorticity (PV)
snapshots after 3 years of integration at low-resolution with equations (1,2)
with (top) standard hyper-viscosity and (bottom) 10 times smaller hyper-
viscosity. We can see aliasing effects on potential vorticity snapshots
integrated with low hyper-viscosity.

###  3.2 Modified viscosity

Here we propose a simple affine modification parameterization of
hyperviscosity. We add a bias to the term $\Delta\mathbf{p}$ in eq. (1), which
becomes $\Delta\left(\mathbf{p}-\mathbf{p}^{\prime}\right)$ where
$\mathbf{p}^{\prime}$ is a dimensional field that does not depend upon time.
The PV advection equation with hyperviscosity becomes

$\displaystyle\partial_{t}\mathbf{q}=\frac{1}{f_{0}}J(\mathbf{q},\mathbf{p})+f_{0}B\mathbf{e}+\frac{1}{f_{0}}\big{(}a_{2}\Delta-
a_{4}\Delta^{2}\big{)}\left(\Delta(\mathbf{p}-\mathbf{p}^{\prime})\right).$
(3)

The elliptic equation (2) remains unchanged.

The goal of this additional term is to reproduce a relevant time-average
pressure field relying on observations or high-resolution solutions. For
example the high-resolution average $\overline{\mathbf{p}_{\rm HR}}$ can be
downsampled to the targeted coarse grid resolution in
$\overline{\mathbf{p}_{\rm HR}}\downarrow$, and we want the average of the
modified low-resolution $\overline{\mathbf{p}_{\rm LR}}$ model to be as close
as possible to the high-resolution reference $\overline{\mathbf{p}_{\rm
HR}}\downarrow$.

We face here an optimal control problem, as the low-resolution average is a
function of the control parameter $\mathbf{p}^{\prime}$. We state it with the
following least-square formulation

$\displaystyle\mathbf{p}^{\prime}_{\rm opt}$
$\displaystyle=\underset{\mathbf{p}^{\prime}}{\text{argmin
}}\mathcal{F}\left(\mathbf{p}^{\prime}\right)$ (4)
$\displaystyle\mathcal{F}\left(\mathbf{p}^{\prime}\right)$
$\displaystyle=\big{\|}\ \overline{\mathbf{p}_{\rm
LR}}\left(\mathbf{p}^{\prime}\right)-\overline{\mathbf{p}_{\rm
HR}}\downarrow\big{\|}^{2}$ (5)

This optimization problem is a priori non-convex and we shall not expect to
find a global optimum. In the following, we propose a numerical procedure to
find a heuristic $\hat{\mathbf{p}^{\prime}}$ of the optimal solution
$\mathbf{p}^{\prime}_{\rm opt}$.

Computationally, the implementation of this modified hyperviscosity is simple
and computationally cheap. We precompute $\Delta\mathbf{p}^{\prime}$ and
subtract it from $\Delta\mathbf{p}$ at each time-integration step. It
increases the integration time of the advection equation (1) by less than 1%
on CPUs and GPUs.

###  3.3 Modified viscosity regularization

The continuously stratified QG equations can be rewritten in a variational
formulation [9] with a Hamiltonian $\mathcal{J}$ defined as

$\mathcal{J}(\mathbf{p})=\frac{1}{2}\int_{\Omega}\frac{1}{f_{0}}|\nabla\mathbf{p}|^{2}+\frac{f_{0}}{N^{2}}(\partial_{z}\mathbf{p})^{2}.$

Our model is a discretized version of the continuous stratification. Since we
add an external wind forcing term and we use an energy conservative Arakawa
advection scheme, we need to add some viscosity or hyperviscosity to dissipate
energy. In a variational formulation, these (hyper-)viscous terms become the
following penalization

$\frac{1}{2}\int_{\Omega}a_{2}|\Delta\mathbf{p}|^{2}+a_{4}\left|\nabla\left(\Delta\mathbf{p}\right)\right|^{2},$

added to the Hamiltonian $\mathcal{J}(\mathbf{p})$ to produce a smooth
solution. The Gradient norm penalization of Laplacian $\mathbf{p}$ guides the
minimization toward solutions of smooth Laplacian. Hyperviscosity corresponds
to the Laplacian norm penalization and enforces a solution of minimum
Laplacian norm. The parameters $a_{2}$ and $a_{4}$ quantify the strength of
these regularization constraints.

Here, we simply propose to replace it with the following penalization

$\frac{1}{2}\int_{\Omega}a_{2}|\Delta(\mathbf{p}-\mathbf{p}^{\prime})|^{2}+a_{4}\left|\nabla\left(\Delta(\mathbf{p}-\mathbf{p}^{\prime})\right)\right|^{2}\
.$

We now penalize $(\mathbf{p}-\mathbf{p}^{\prime})$ instead of $\mathbf{p}$,
meaning that we guide the solution to a possibly non-smooth reference
$\mathbf{p}^{\prime}$ that will produce the correct large scale behavior.

###  3.4 Iterative procedure

Here we present a method to find a solution to the optimization problem (4). A
natural guess for $\mathbf{p}^{\prime}_{\rm opt}$ is
$\overline{\mathbf{p}_{\rm HR}}\downarrow$. We solve the equations and compute
the average pressure $\overline{\mathbf{p}_{\rm LR}}$. Results are shown in
Fig. 4. It is a good first-guess, but the difference
$\overline{\mathbf{p}_{\rm HR}}\downarrow-\overline{\mathbf{p}_{\rm LR}}$ is
still large.

We propose the following iterative procedure to find a better guess for
$\mathbf{p}^{\prime}_{\rm opt}$. In the following we assume that we are in low
resolution, i.e. $\mathbf{p}=\mathbf{p}_{\rm LR}$ and
$\overline{\mathbf{p}}=\overline{\mathbf{p}_{\rm LR}}$ unless explicitly
written.

  * •

We set $\mathbf{p}^{\prime}_{0}$ and we compute the average pressure
$\overline{\mathbf{p}}_{0}$ solving standard equations (1,2) without modified
viscosity.

  * •

Choose $k\in\left]0,1\right]$.

  * •

Start with $\mathbf{p}^{\prime}_{1}=\overline{\mathbf{p}_{\rm HR}}\downarrow$.

  * •

Evolve the ensemble for $n$ years and compute the corresponding average
pressure $\overline{\mathbf{p}}_{1}$ with ensemble average.

  * •

For $n=1\dots$:

    * –

Set
$\mathbf{p}^{\prime}_{n+1}=\mathbf{p}^{\prime}_{n}+k\left(\overline{\mathbf{p}_{\rm
HR}}\downarrow-\overline{\mathbf{p}}_{n}\right)$.

    * –

Evolve the ensemble for $n$ years and compute new average pressure
$\overline{\mathbf{p}}_{n+1}$.

  * •

return $\mathbf{p}^{\prime}_{n}$ and $\overline{\mathbf{p}}_{n}$

There is no theoretical guarantee that this procedure converges, but we
observe in the next section that it converges with the double-gyre QG model
that we use.

##  4 Results and discussion

###  4.1 Statistics

We use ensemble averages to compute the statistics. To create ensembles of
size $N$, we start from a zero solution and spin up the models for 100 years
with a timestep of 1200 seconds to reach statistically steady states as in
[14]. Then we run the models for 500 years and save 10 snapshots a year to get
5000 snapshots, and we randomly select $N$ snapshots out of these 5000
snapshots. The ensemble averages are simply average over these $N$ ensemble
members that we evolve in parallel. Such ensemble averages are denoted with
$\overline{\bullet}$ in the following, i.e. the average pressure is denoted by
$\overline{\mathbf{p}}$, average velocity by $\overline{u}$, etc.

###  4.2 Iterative procedure

Figure 3: Evolution of the relative square error
$\frac{\|\overline{\mathbf{p}}_{n}-\overline{\mathbf{p}_{HR}}\downarrow\|^{2}}{\|\overline{\mathbf{p}_{HR}}\downarrow\|^{2}}$
w.r.t iterations of the procedure.

We test the iterative procedure described in section 4.2 with the double-gyre
model presented in section 2 in the eddy-permitting regime. We use $n=10$
years to evolve the ensemble after each iterate. We compute the reference
pressure average $\overline{\mathbf{p}}_{\rm HR}$ with the same model in the
eddy-resolving regime.

Figure 4: Top-layer average pressure (top) and velocity (bottom) of (left-to-
right) proposed model at low-resolution, reference, and the difference between
the two.

Fig. 3 shows the relative square error
$\|\overline{\mathbf{p}}_{n}-\overline{\mathbf{p}_{HR}}\downarrow\|^{2}/\|\overline{\mathbf{p}_{HR}}\downarrow\|^{2}$
at iterations of the procedure with $k=1$ and $k=0.7$. The procedure converges
with $k=0.7$ and oscillates with $k=1$.

Figure 5: Top-layer kinetic-energy spectra average with models at high-
resolution (HR), at low-resolution (LR) and at low-resolution with proposed
modified viscosity. The decreasing slope of the spectrum of the proposed model
is much closer to the high-resolution reference.

Fig. 4 shows the output average pressure $\overline{\mathbf{p}}_{n}$ of the
iterative procedure, the reference $\overline{\mathbf{p}}_{\rm HR}$ and the
difference between the two, as well as for zonal velocity $\mathbf{u}$ . Our
model can reproduce the eastward jet produced by the high-resolution reference
model. Kinetic energy spectra shown on Fig. 5 shows also the improvement of
our model compared to low-resolution. Finally, Fig. 6 shows high-resolution
and low-resolution snapshots as well as a snapshot of the proposed model at
low-resolution. Our model effectively produces the eastward jet and a re-
circulation zone around it where eddies are created. Artifacts can be also
observed on the zonal velocity and potential vorticity on the right of Fig. 6.
They can likely be Rossby waves created by the harmonic regularization terms,
which remain an artificial constraint, but this needs to be studied further.

Figure 6: PV and zonal velocity snapshots form (left-to-right) high-
resolution, low-resolution and proposed model at low-resolution.

##  5 Conclusion

We presented a simple modified-viscosity scheme for coarse resolution ocean
modeling that we derived and tested on a double-gyre multi-layer quasi-
geostrophic model. We interpret it as a modified regularization technique that
will guide the solution to a reference rather than producing a too smooth
solution in the eddy-permitting regime. The technique requires solving an
optimization problem, and we presented a procedure to find a good guess for
the solutions. We showed that it converges to a reasonable solution that
fairly reproduces the input reference.

If this method mimics the average of the high-resolution, it only reproduces
partially the variability and higher-order statistics of the high-resolution.
We see in Fig. 5 our model’s snapshots resemble the averages. In future works,
we consider using this method as a deterministic basis for stochastic
parameterizations such as Location-Uncertainty [12].

## Acknowledgments

The authors acknowledge the support of the ERC EU project 856408-STUOD.

## Appendix

### Downsampling procedure

Downsampling the high-resolution solution on a low-resolution grid consists of
interpolating the high-resolution (769$\times$ 961) streamfunction on the low-
resolution (97$\times$121) grid. Then we can compute the potential vorticity
using equation (2). Because of the no-flow constraint, the downsampled
streamfunction should be constant on the boundaries and should satisfy a mass
conservation constraint [8]. We also want to preserve the frequency
information and prevent aliasing.

We use the following procedure:

  1. 1.

we apply a Gaussian filter and downsample the streamfunction on the domain
except on the boundaries.

  2. 2.

we adding homogeneous solutions of the Elliptic equation (2) to the
streamfunction in order to satisfy the mass conservation as in [8].

### Parameter tables

Parameters | Value | Description
---|---|---
$L_{x}\times L_{y}$ |  $(3840\times 4800)$ km | Domain size
$H_{k}$ |  $(350,750,2900)$ m | Mean layer thickness
$g_{k}^{\prime}$ |  $(0.025,0.0125)$ m s-2 | Reduced grativity
$\delta_{\text{ek}}$ |  $2$ m | Bottom Ekman layer thickness
$\tau_{0}$ |  $2\times 10^{-5}$ m2 s-2 | Wind stress magantitude
$a_{2}$ |  $0$ m2 s-1 | Laplacian viscosity coefficient
$f_{0}$ |  $9.375\times 10^{-5}$ s-1 | Mean Coriolis parameter
$\beta$ |  $1.754\times 10^{-11}$ (m s)-1 | Coriolis parameter gradient
Table 1: Common parameters for all the models. Grid dimensions | Resolution | Timestep | Hyperviscosity ($a_{4}$)
---|---|---|---
769 $\times$ 961 |  $5$ km |  $600$ s |  $2.0\times 10^{9}$ m4 s-1
97 $\times$ 121 |  $40$ km |  $1200$ s |  $5.0\times 10^{11}$ m4 s-1
Table 2: Grid-dependent parameters.

## References

  * [1] Arakawa, Akio and Lamb, Vivian R (1981), A potential enstrophy and energy conserving scheme for the shallow water equations. Monthly Weather Review 109:18–36.
  * [2] Fox-Kemper, B and Menemenlis, D (2008), Can large eddy simulation techniques improve mesoscale rich ocean models?. Washington DC American Geophysical Union Geophysical Monograph Series 177:319–337.
  * [3] Fox-Kemper, B and Bachman, S and Pearson, B and Reckinger, S (2014), Principles and advances in subgrid modelling for eddy-rich simulations. CLIVAR Exchanges No. 65, Vol. 19, No. 2.
  * [4] Frisch, Uriel and Kurien, Susan and Pandit, Rahul and Pauls, Walter and Ray, Samriddhi Sankar and Wirth, Achim and Zhu, Jian-Zhou (2008), Hyperviscosity, Galerkin truncation, and bottlenecks in turbulence. Physical review letters 101:144501.
  * [5] Griffies, Stephen M and Hallberg, Robert W (2000), Biharmonic friction with a Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models. Monthly Weather Review 128:2935–2946.
  * [6] Hogg, Andrew Mc C and Dewar, William K and Killworth, Peter D and Blundell, Jeffrey R (2003), A quasi-geostrophic coupled model (Q-GCM). Monthly Weather Review 131:2261–2278.
  * [7] Hogg, Andrew Mc C and Killworth, Peter D and Blundell, Jeffrey R and Dewar, William K (2005), Mechanisms of decadal variability of the wind-driven ocean circulation. Journal of physical oceanography 35:512–531.
  * [8] Hogg, AM and Blundell, JR and Dewar, WK and Killworth, PD (2014), Formulation and users’ guide for Q-GCM, http://q-gcm.org/downloads/q-gcm-v1.5.0.pdf.
  * [9] Darryl D. Holm, Vladimir Zeitlin (1998), Hamilton’s principle for quasigeostrophic motion. Physics of fluids 10(4):800–806.
  * [10] Jansen, Malte F and Held, Isaac M (2014), Parameterizing subgrid-scale eddy effects using energetically consistent backscatter. Ocean Modelling 80:36–48.
  * [11] Leith, CE (1996), Stochastic models of chaotic systems. Physica D: Nonlinear Phenomena 98:481–491.
  * [12] Mémin, Etienne (2014), Fluid flow dynamics under location uncertainty. Geophysical & Astrophysical Fluid Dynamics 108:119–146.
  * [13] Paszke, Adam and Gross, Sam and Massa, Francisco and Lerer, Adam and Bradbury, James and Chanan, Gregory and Killeen, Trevor and Lin, Zeming and Gimelshein, Natalia and Antiga, Luca and others (2019), Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems 32:8026–8037.
  * [14] Ryzhov, EA and Kondrashov, D and Agarwal, N and McWilliams, JC and Berloff, P (2020), On data-driven induction of the low-frequency variability in a coarse-resolution ocean model. Ocean Modelling 153:101664
  * [15] Roullet, Guillaume and Gaillard, Tugdual (2021), A Fast Monotone Discretization of the Rotating Shallow Water Equations, ESSOAR/JAMES.
  * [16] Shevchenko, IV and Berloff, PS (2015), Multi-layer quasi-geostrophic ocean dynamics in eddy-resolving regimes. Ocean Modelling 94:1–14.
  * [17] Smagorinsky J. (1963), General circulation experiments with the primitive equations. Monthly weather review 91.3:99-164.
  * [18] Vallis, G. K. (2017), Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation. Cambridge University Press.


#  $R$-matrix for the XX spin chain at the special imaginary value of
staggered magnetic field

P. N. Bibikov

(Russian State Hydrometeorological University, Saint-Petersburg, Russia)

###### Abstract

It is pointed, that the $16\times 16$ Hamiltonian density matrix,
corresponding to XX spin chain in a staggered magnetic field, satisfies the
Reshetikhin condition as well, as the two higher ones. All of them are
necessary for the existence of the corresponding $R$-matrix. At the special
imaginary value of the staggered field, the $R$-matrix is obtained explicitly.
It is shown that the corresponding Hamiltonian has purely imaginary spectrum.

The XX model, suggested in [1], and related to the Hamiltonian

$\hat{H}=\sum_{n}H_{n,n+1},$ (1)

where

$H_{n,n+1}=\frac{1}{2}\Big{(}{\bf S}^{+}_{n}{\bf S}^{-}_{n+1}+{\bf
S}^{-}_{n}{\bf S}^{+}_{n+1}\Big{)},$ (2)

is now one of the most studied integrable spin chains (see [2] and references
therein). Here ${\bf S}^{j}_{n}$ ($j=+,-,z$) are the spin-1/2 operators
associated with the $n$-th site of the chain.

The integrability of (1), (2) has been first demonstrated within the Jordan-
Wigner transformation to the free-fermion model and then confirmed by various
alternative approaches (see [2] and references therein). Namely, the
Hamiltonian density $4\times 4$ matrix

$H=\frac{1}{2}\Big{(}{\bf S}^{+}\otimes{\bf S}^{-}+{\bf S}^{-}\otimes{\bf
S}^{+}\Big{)},$ (3)

related to (2), may be represented in the form

$H=\frac{1}{2}\frac{dR(\lambda)}{d\lambda}\Big{|}_{\lambda=0},$ (4)

where the $R$-matrix

$R(\lambda)=\left(\begin{array}[]{cccc}\cos{\lambda}&0&0&0\\\
0&1&\sin{\lambda}&0\\\ 0&\sin{\lambda}&1&0\\\
0&0&0&\cos{\lambda}\end{array}\right),$ (5)

satisfies the Yang-Baxter equation (in the Braid group form) [4]

$R_{12}(\lambda-\mu)R_{23}(\lambda)R_{12}(\mu)=R_{23}(\mu)R_{12}(\lambda)R_{23}(\lambda-\mu),$
(6)

(as usual, $R_{12}(\lambda)\equiv R(\lambda)\otimes I_{2}$ and
$R_{23}(\lambda)\equiv I_{2}\otimes R(\lambda)$) and the initial condition

$R(0)=I_{4}.$ (7)

Here by $I_{M}$ we denote the $M\times M$ identity matrix.

Within the slightly modified version of the [1] approach, it has been proved
[3], that the model (1), (2) remains integrable even under addition of a
staggered longitudinal magnetic field, when

$H_{n,n+1}(h_{\rm st})=\frac{1}{2}\Big{[}{\bf S}^{+}_{n}{\bf S}^{-}_{n+1}+{\bf
S}^{-}_{n}{\bf S}^{+}_{n+1}+h_{\rm st}\Big{(}(-1)^{n}{\bf
S}^{z}_{n}+(-1)^{n+1}{\bf S}^{z}_{n+1}\Big{)}\Big{]},$ (8)

instead of (2). Namely, with the use of the Fermi operators representation

${\bf S}^{+}_{n}=\prod_{m<n}\sigma_{m}a_{n},\qquad{\bf
S}^{-}_{n}=\prod_{m<n}\sigma_{m}a_{n}^{\dagger},\qquad{\bf
S}^{z}_{n}=\frac{1}{2}[a_{n},a_{n}^{\dagger}]=\frac{1}{2}-a_{n}^{\dagger}a_{n},$
(9)

where $\sigma_{n}\equiv 1-2a_{n}^{\dagger}a_{n}$ and $\\{a_{m},a_{n}\\}=0$,
and $\\{a_{m},a_{n}^{\dagger}\\}=\delta_{mn}$, the Hamiltonian (1), (2) takes
the well tractable bilinear form

$\hat{H}(h_{\rm
st})=\frac{1}{2}\sum_{n}\Big{[}a_{n+1}^{\dagger}a_{n}+a_{n}^{\dagger}a_{n+1}-h_{\rm
st}\Big{(}(-1)^{n}a_{n}^{\dagger}a_{n}+(-1)^{n+1}a_{n+1}^{\dagger}a_{n+1}\Big{)}\Big{]},$
(10)

which, under the substitutions $a_{2n}\equiv b_{n}$, and
$a_{2n+1}\equiv\tilde{b}_{n}$, turns into

$\hat{H}(h_{\rm
st})=\frac{1}{2}\sum_{n}\Big{(}\tilde{b}_{n}^{\dagger}b_{n}+b_{n}^{\dagger}\tilde{b}_{n}+b_{n+1}^{\dagger}\tilde{b}_{n}+\tilde{b}_{n}^{\dagger}b_{n+1}+h_{\rm
st}(2\tilde{b}_{n}^{\dagger}\tilde{b}_{n}-b_{n}^{\dagger}b_{n}-b_{n+1}^{\dagger}b_{n+1})\Big{)}.$
(11)

With the use of the Fourier representations

$b_{n}=\frac{1}{\sqrt{2\pi}}\int_{0}^{2\pi}dk{\rm
e}^{-ikn}c_{1}(k),\qquad\tilde{b}_{n}=\frac{1}{\sqrt{2\pi}}\int_{0}^{2\pi}dk{\rm
e}^{-ikn}c_{2}(k),$ (12)

where

$\\{c_{j}^{\dagger}(k),c_{l}(p)\\}=\delta_{jl}\delta(k-p),\qquad\\{c_{j}(k),c_{l}(p)\\}=0,$
(13)

(10) is reduced to the form

$\hat{H}(h_{\rm
st})=\int_{0}^{2\pi}dk\left(\begin{array}[]{cc}c_{1}^{\dagger}(k)&c_{2}^{\dagger}(k)\end{array}\right)M(k,h_{\rm
st})\left(\begin{array}[]{c}c_{1}(k)\\\ c_{2}(k)\end{array}\right),$ (14)

where

$M(k,h_{\rm st})=\frac{1}{2}\left(\begin{array}[]{cc}-2h_{\rm st}&1+{\rm
e}^{ik}\\\ 1+{\rm e}^{-ik}&2h_{\rm st}\end{array}\right).$ (15)

As it has been shown in [3], for real $h_{\rm st}$ the Hamiltonian (14) may be
diagonalized by unitary transformation, so that the model (1), (2) is
integrable.

It is desirable to confirm this result by presentation of the corresponding
$R$-matrix. Since (8) is invariant only under even translations
($H_{n,n+1}=H_{n+2,n+3}\neq H_{n+1,n+2}$), the approach based on a $4\times 4$
matrix needs a modification. Namely, one should pass from (3) to the $16\times
16$ matrix, corresponding to the representation

$\hat{H}(h_{\rm st})=\sum_{n}{\cal H}_{2n,2n+3}(h_{\rm st}),$ (16)

more general than (1). Here

${\cal H}_{2n,2n+3}(h_{\rm st})=\frac{1}{2}\Big{(}H_{2n,2n+1}(h_{\rm
st})+H_{2n+2,2n+3}(h_{\rm st})\Big{)}+H_{2n+1,2n+2}(h_{\rm st}),$ (17)

and the corresponding $16\times 16$ Hamiltonian density matrix is

${\cal H}(h_{\rm st})=\frac{1}{2}\Big{(}H^{(+)}(h_{\rm st})\otimes
I_{4}+I_{4}\otimes H^{(+)}(h_{\rm st})\Big{)}+I_{2}\otimes H^{(-)}(h_{\rm
st})\otimes I_{2},$ (18)

where

$H^{(\pm)}(h_{\rm st})=\frac{1}{2}\Big{[}{\bf S}^{+}\otimes{\bf S}^{-}+{\bf
S}^{-}\otimes{\bf S}^{+}\pm h_{\rm st}\Big{(}{\bf S}^{z}\otimes
I_{2}-I_{2}\otimes{\bf S}^{z}\Big{)}\Big{]}.$ (19)

The direct calculation shows, that the matrix (18) satisfies the Reshetikhin
condition [4]

$[{\cal H}(h_{\rm st})\otimes I_{4}+I_{4}\otimes{\cal H}(h_{\rm st}),[{\cal
H}(h_{\rm st})\otimes I_{4},I_{4}\otimes{\cal H}(h_{\rm
st})]]=I_{4}\otimes{\cal K}-{\cal K}\otimes I_{4},$ (20)

(where the $16\times 16$ matrix ${\cal K}$ may be readily calculated).

According to (20), the sum

$\hat{Q}=\sum_{n}[{\cal H}_{2n-2,2n+1}(h_{\rm st}),{\cal H}_{2n,2n+3}(h_{\rm
st})],$ (21)

is the first integral of motion

$[\hat{H},\hat{Q}]=0.$ (22)

In fact, $\hat{Q}$ does not depend on $h_{\rm st}$ and may be reduced to the
form

$\hat{Q}=\frac{1}{2}\sum_{n}[H_{n-1,n},H_{n,n+1}].$ (23)

Really, following (17),

$[{\cal H}_{2n-2,2n+1}(h_{\rm st}),{\cal H}_{2n,2n+3}(h_{\rm
st})]=\frac{1}{2}\Big{(}[H_{2n-1,2n}^{(-)},H^{(+)}_{2n,2n+1}]+[H^{(+)}_{2n,2n+1},H^{(-)}_{2n+1,2n+2}]\Big{)}.$
(24)

At the same time,

$\displaystyle[H^{(+)}_{2n,2n+1},H^{(-)}_{2n+1,2n+2}]=[H_{2n,2n+1},H_{2n+1,2n+2}]+h_{\rm
st}[{\bf S}^{z}_{2n+1},H_{2n,2n+1}-H_{2n+1,2n+2}],$
$\displaystyle[H_{2n-1,2n}^{(-)},H_{2n,2n+1}^{(+)}]=[H_{2n-1,2n},H_{2n,2n+1}]+h_{\rm
st}[{\bf S}^{z}_{2n},H_{2n,2n+1}-H_{2n-1,2n}].$ (25)

Since $[{\bf S}^{z}_{n}+{\bf S}^{z}_{n+1},H_{n,n+1}]=0$, the substitution of
($R$-matrix for the XX spin chain at the special imaginary value of staggered
magnetic field) into (24) gives

$\displaystyle[{\cal H}_{2n-2,2n+1}(h_{\rm st}),{\cal H}_{2n,2n+3}(h_{\rm
st})]=\frac{1}{2}\Big{(}[H_{2n-1,2n},H_{2n,2n+1}]+[H_{2n,2n+1},H_{2n+1,2n+2}]\Big{)}$
$\displaystyle+\frac{h_{\rm st}}{4}\Big{(}[{\bf S}^{z}_{2n-1}-{\bf
S}^{z}_{2n},H_{2n-1,2n}]-[{\bf S}^{z}_{2n+1}-{\bf
S}^{z}_{2n+2},H_{2n+1,2n+2}]\Big{)}.$ (26)

Under the substitution of ($R$-matrix for the XX spin chain at the special
imaginary value of staggered magnetic field) into (21), the terms proportional
to $h_{\rm st}$ are canceled, and one gets (23).

It may be readily proved by direct calculations, that additionally to (20) the
Hamiltonian density matrix satisfies two higher order integrability
conditions, presented in [5] and [6]. Both of them also are necessary (but not
sufficient) for representation of ${\cal H}(h_{\rm st})$ in the form (4), but
with the $16\times 16$ $R$-matrix ${\cal R}(\lambda,h_{\rm st})$, which
satisfies as the Yang-Baxter equation (6), as the initial condition (7) (with
$I_{4}$ replaced by $I_{16}$).

As it has been suggested in [6], and shown in [7, 8], in such a situation one
may try to guess the entries of $R$-matrix by analyzing the series expansions
of their ratios. For the model (16) the author succeed, however, only in the
two special cases, $h_{\rm st}=0$ and $h_{\rm st}=i$.

The direct analysis of power series for the entries shows, that both the
$R$-matrices have the similar structure

${\cal
R}=\left(\begin{array}[]{cccccccccccccccc}*&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&*&*&0&*&0&0&0&0&0&0&0&0&0&0&0\\\ 0&*&*&0&*&0&0&0&*&0&0&0&0&0&0&0\\\
0&0&0&*&0&*&*&0&0&*&*&0&*&0&0&0\\\ 0&*&*&0&*&0&0&0&*&0&0&0&0&0&0&0\\\
0&0&0&*&0&*&*&0&0&*&*&0&*&0&0&0\\\ 0&0&0&*&0&*&*&0&0&*&*&0&*&0&0&0\\\
0&0&0&0&0&0&0&*&0&0&0&*&0&*&0&0\\\ 0&0&*&0&*&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&*&0&*&*&0&0&*&*&0&*&0&0&0\\\ 0&0&0&*&0&*&*&0&0&*&*&0&*&0&0&0\\\
0&0&0&0&0&0&0&*&0&0&0&*&0&*&*&0\\\ 0&0&0&*&0&*&*&0&0&*&*&0&*&0&0&0\\\
0&0&0&0&0&0&0&*&0&0&0&*&0&*&*&0\\\ 0&0&0&0&0&0&0&0&0&0&0&*&0&*&*&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&*\\\ \end{array}\right),$ (27)

and for both of them

${\cal R}_{ji}={\cal R}_{ij},\qquad{\cal R}_{17-i,17-j}=\bar{\cal R}_{ij},$
(28)

(for the shortness we put ${\cal R}_{ij}\equiv{\cal R}_{ij}(\lambda,h_{\rm
st})$). Moreover,

$\displaystyle{\cal R}_{16,16}={\cal R}_{11},\quad{\cal R}_{15,15}={\cal
R}_{99},\quad{\cal R}_{14,14}={\cal R}_{55},$ $\displaystyle{\cal
R}_{12,12}={\cal R}_{33},\quad{\cal R}_{10,10}={\cal R}_{77},\quad{\cal
R}_{8,8}={\cal R}_{22},$ $\displaystyle{\cal R}_{14,15}={\cal
R}_{59},\quad{\cal R}_{10,11}={\cal R}_{7,11},\quad{\cal R}_{10,13}={\cal
R}_{7,13},$ $\displaystyle{\cal R}_{8,12}={\cal R}_{23},\quad{\cal
R}_{6,10}={\cal R}_{67},\quad{\cal R}_{4,10}={\cal R}_{47},$
$\displaystyle{\cal R}_{7,10}={\cal R}_{39}={\cal R}_{8,14}={\cal
R}_{12,15}={\cal R}_{25}={\cal R}_{6,11}-{\cal R}_{4,13}={\cal R}_{77}-{\cal
R}_{11},$ $\displaystyle{\cal R}_{12,14}={\cal R}_{35}={\cal R}_{67}+{\cal
R}_{10,11}.$ (29)

According to (28), it convenient to introduce the real variables

${\cal A}_{ij}=\frac{{\cal R}_{ij}+{\cal R}_{17-i,17-j}}{2},\qquad{\cal
B}_{ij}=\frac{{\cal R}_{ij}-{\cal R}_{17-i,17-j}}{2i},$ (30)

for which

${\cal A}_{ji}={\cal A}_{ij},\quad{\cal B}_{ji}={\cal B}_{ij},\quad{\cal
A}_{17-i,17-j}={\cal A}_{ij},\quad{\cal B}_{17-i,17-j}=-{\cal B}_{ij}.$ (31)

Following ($R$-matrix for the XX spin chain at the special imaginary value of
staggered magnetic field),

$\displaystyle{\cal A}_{55}={\cal A}_{33},\quad{\cal A}_{88}={\cal
A}_{22},\quad A_{4,10}=A_{47},\quad{\cal A}_{59}={\cal A}_{23},\quad{\cal
A}_{6,10}={\cal A}_{67},$ $\displaystyle{\cal A}_{39}={\cal A}_{25}={\cal
A}_{7,10}={\cal A}_{77}-{\cal A}_{11},\quad{\cal A}_{4,13}={\cal
A}_{6,11}+{\cal A}_{11}-{\cal A}_{77},$ $\displaystyle{\cal A}_{7,11}={\cal
A}_{35}-A_{67},$ (32)

and

$\displaystyle{\cal B}_{11}={\cal B}_{77}={\cal B}_{25}={\cal B}_{35}={\cal
B}_{39}=0,\quad{\cal B}_{55}=-{\cal B}_{33},\quad{\cal B}_{88}={\cal B}_{22},$
$\displaystyle{\cal B}_{59}=-{\cal B}_{23},\quad{\cal B}_{4,10}={\cal
B}_{47},\quad{\cal B}_{6,10}={\cal B}_{67},$ (33)

According to (30)-($R$-matrix for the XX spin chain at the special imaginary
value of staggered magnetic field), all the entries of the $R$-matrix ${\cal
R}(\lambda,h_{\rm st})$ (27) are linear combinations of 13 parameters
$\\{{\cal A}_{11},{\cal A}_{22},{\cal A}_{33},{\cal A}_{44},{\cal
A}_{66},{\cal A}_{77},{\cal A}_{23},{\cal A}_{35},{\cal A}_{46},{\cal
A}_{47},{\cal A}_{4,11},{\cal A}_{67},{\cal A}_{6,11}\\}$, and 9 parameters
$\\{{\cal B}_{22},{\cal B}_{33},{\cal B}_{44},{\cal B}_{66},{\cal
B}_{23},{\cal B}_{46},{\cal B}_{47},{\cal B}_{4,11},{\cal B}_{67}\\}$. Namely,
for $h_{\rm st}=0$ one may put

$\displaystyle{\cal A}_{11}={\tt c}^{4},\quad{\cal A}_{22}={\tt
c}^{3},\quad{\cal A}_{33}={\tt c}(1+{\tt s}^{2}),\quad{\cal
A}_{44}=1,\quad{\cal A}_{66}=1+{\tt s}^{2}-{\tt s}^{4},$ $\displaystyle{\cal
A}_{77}={\tt c}^{2},\quad{\cal A}_{23}={\tt c}^{2}{\tt s},\quad{\cal
A}_{35}=2{\tt c}{\tt s},\quad{\cal A}_{46}=(2-{\tt s}^{2}){\tt s},\quad{\cal
A}_{47}={\tt c}{\tt s}^{2},\quad{\cal A}_{4,11}={\tt s}^{3},$
$\displaystyle{\cal A}_{67}={\tt c}{\tt s},\quad{\cal A}_{6,11}={\tt
s}^{2},\quad{\cal B}_{ij}=0,\qquad{\tt c}\equiv\cos{\lambda},\quad{\tt
s}\equiv\sin{\lambda}.$ (34)

At the same time, for $h_{\rm st}=i$

$\displaystyle{\cal A}_{11}={\cal A}_{22}={\tt ch}^{4},\quad{\cal
A}_{33}=(2-{\tt ch}^{2}){\tt ch}^{2},\quad{\cal A}_{44}=2{\tt
ch}^{2}-1,\quad{\cal A}_{66}=1-3{\tt ch}^{2}{\tt sh}^{2},$ $\displaystyle{\cal
A}_{77}=(2{\tt ch}^{2}-1){\tt ch}^{2},\quad{\cal A}_{23}={\tt ch}^{3}{\tt
sh},\quad{\cal A}_{35}=2{\tt ch}{\tt sh},\quad{\cal A}_{46}=(3-{\tt
ch}^{2}){\tt ch}{\tt sh},$ $\displaystyle{\cal A}_{47}={\tt ch}^{2}{\tt
sh}^{2},\quad{\cal A}_{4,11}={\tt ch}{\tt sh}^{3},\quad{\cal A}_{67}={\tt
ch}{\tt sh},\quad{\cal A}_{6,11}=(2{\tt ch}^{2}-1){\tt sh}^{2},$
$\displaystyle{\cal B}_{22}={\tt ch}^{3}{\tt sh},\quad{\cal B}_{33}=-(2+{\tt
ch}^{2}){\tt ch}{\tt sh},\quad{\cal B}_{44}=-2{\tt ch}{\tt sh},\quad{\cal
B}_{66}=4{\tt ch}{\tt sh},$ $\displaystyle{\cal B}_{23}=-{\tt ch}^{2}{\tt
sh}^{2},\quad{\cal B}_{46}=(1+{\tt ch}^{2}){\tt sh}^{2},\quad{\cal
B}_{47}={\tt ch}{\tt sh}^{3},\quad{\cal B}_{4,11}=-{\tt sh}^{4}$
$\displaystyle{\cal B}_{67}=2{\tt ch}^{2}{\tt sh}^{2},\qquad{\tt
ch}\equiv\cosh{\lambda},\quad{\tt sh}\equiv\sinh{\lambda}.$ (35)

In both the cases the correspondence between ${\cal H}$ and ${\cal R}$ is
given exactly by the formula (4) (with $H,R\rightarrow{\cal H},{\cal R}$).

Let us study the case $h_{\rm st}=i$ in more details. Introducing the
operators

$\displaystyle c_{+}(k)=\frac{c_{1}(k)+\gamma
c_{2}(k)}{\sqrt{1+|\gamma|^{2}}},\qquad
c_{-}(k)=\frac{\bar{\gamma}c_{1}(k)+c_{2}(k)}{\sqrt{1+|\gamma|^{2}}},\qquad\gamma\equiv\frac{1+{\rm
e}^{ik}}{2i(1-\sin{k/2})},$ (36)

and using the identities

$|\gamma|^{2}=\frac{\cos^{2}{k/2}}{(1-\sin{k/2})^{2}},\qquad
1+|\gamma|^{2}=\frac{2}{1-\sin{k/2}},\qquad
1-|\gamma|^{2}=-\frac{2\sin{k/2}}{1-\sin{k/2}},$ (37)

one readily gets from (13) and (14)

$\\{c_{\pm}^{\dagger}(k),c_{\pm}(p)\\}=\delta(k-p),\qquad\\{c_{+}^{\dagger}(k),c_{-}(p)\\}=i\frac{1+{\rm
e}^{-ik}}{2}\delta(k-p),$ (38)

and

$\displaystyle\hat{H}(i)=\int_{0}^{2\pi}\frac{dk}{\sin{k/2}}\Big{[}i(c_{+}^{\dagger}(k)c_{+}(k)-c_{-}^{\dagger}(k)c_{-}(k))$
$\displaystyle-\frac{1}{2}\Big{(}(1+{\rm
e}^{ik})c_{+}^{\dagger}(k)c_{-}(k)+(1+{\rm
e}^{-ik})c_{-}^{\dagger}(k)c_{+}(k)\Big{)}\Big{]}.$ (39)

According to (38) and ($R$-matrix for the XX spin chain at the special
imaginary value of staggered magnetic field),

$[\hat{H}(i),c^{\dagger}_{\pm}(k)]=\pm
i\sin{\frac{k}{2}}c^{\dagger}_{\pm}(k),$ (40)

and

$\displaystyle[\hat{H}(i),c_{+}(k)]=\frac{1}{\sin{k/2}}\Big{[}-i\Big{(}1+\cos^{2}{\frac{k}{2}}\Big{)}c_{+}(k)+(1+{\rm
e}^{ik})c_{-}(k)\Big{]},$
$\displaystyle[\hat{H}(i),c_{-}(k)]=\frac{1}{\sin{k/2}}\Big{[}(1+{\rm
e}^{-ik})c_{+}(k)+i\Big{(}1+\cos^{2}{\frac{k}{2}}\Big{)}c_{-}(k)\Big{]}.$ (41)

As it follows from (40) and ($R$-matrix for the XX spin chain at the special
imaginary value of staggered magnetic field), the spectrum of $\hat{H}(i)$
should be purely imaginary. Hence, the spectrum of $i\hat{H}(i)$ should be
purely real.

At the present time an interest to the non-Hermitian quantum spin chains (and
especially with real spectrums) has steadily grown (see [9, 10, 11, 12] and
referenced therein). The author believes, that the suggested model also may be
interesting.

## References

  * [1] Lieb E, Schultz T, and Mattis D 1961 Two soluble models of an antiferromagnetic chain, Ann. Phys. 16, 407-466
  * [2] ${\rm G\ddot{o}hmann}$ F, Kozlowski K K, Sirker J, Suzuki J 2019 Equilibrium dynamics of the XX chain, Phys. Rev. B 100, 155428
  * [3] Kontorovich V M and Tsukernik V M 1968 Magnetic properties of a spin array with two sublattices, Soviet Physics JETP, 26, 687-691
  * [4] Kulish P P, Sklyanin E K 1982 Quantum spectral transform method. Recent developments, Lecture Notes in Physics 151 ed. J. Hietarinta and C. Montonen, Springer (1982)
  * [5] $\rm M\ddot{u}tter$ K H and Shmitt P 1995 Solvable spin-1 models in one dimension, J. Phys. A: Math. Gen. 28 2265-2276
  * [6] Bibikov P N 2003 How to solve Yang-Baxter equation using the Taylor expansion of $R$-matrix, Phys. Lett. A 314 209-213
  * [7] Bibikov P N 2005 $R$-matrices for $SU(2)$-invariant two-leg spin ladders, J. Math. Sci. (N. Y.) 125 111-117
  * [8] Bibikov P N 2007 $R$-matrices for integrable $SU(2)\times U(1)$-symmetric $S=1/2$ spin-orbital chains, J. Phys. A: Math. and Theor. 40 4683
  * [9] $\rm Bu\check{c}a$ B, Booker C, Medenjak M, and Jaksch D 2020 Bethe ansatz approach for dissipation: exact solutions of quantum many-body dynamics under loss, New. Journ. Phys. 22 123040
  * [10] Chen G, Song F, Lado J L 2023 Topological spin excitations in non-Hermitian spin chains with a generalized kernel polynomial algorithm, Phys. Rev. Lett. 130 100401
  * [11] Sarkar A, Sen S, Kumar S 2023 Spectral crossover in non-Hermitian spin chains: Comparison with random matrix theory, Phys. Rev. E 108 054210
  * [12] Muldoon J, Joglekar Y N 2023 ${\cal PT}$-symmetry breaking in quantum spin chains with exceptional non-Hermicities, Phys. Rev. A 108 062205