# Effect of junction critical current disorder in superconducting quantum
interference filter arrays

K.-H. Müller<EMAIL_ADDRESS>E. E. Mitchell  CSIRO Manufacturing, PO Box
218, Lindfield, NSW 2070, Australia

###### Abstract

In this study, we investigated the performance of two 2D superconducting
quantum interference filter (SQIF) arrays fabricated from YBCO thin films at a
temperature of 77 K. Each array consisted of 6 Josephson junctions (JJs) in
parallel and 167 in series. We conducted both experimental and theoretical
analyses, measuring the arrays’ voltage responses to an applied magnetic field
and their voltage versus bias-current characteristics. To properly model the
planar array layouts, our theoretical model used the stream function approach
and also included the Johnson noise in the JJs. The model further divides the
superconducting current density of the arrays into its Meissner current,
circulating current, and bias current parts for practicality. Since the
fabrication process of YBCO thin films cannot produce identical JJ critical
currents, we assumed a log-normal distribution to model the JJ critical
current disorder. Our model predictions, with a JJ critical current spread of
50%, agreed well with our experimental data. Using our model, we were able to
study the dependence of the voltage modulation depth on critical current
disorder and London penetration depth. We also analyzed the observed
reflection asymmetries of the voltage versus magnetic field characteristics,
which might provide insight into the degree of critical current disorder.
Overall, our findings suggest that the use of YBCO thin films in SQIF arrays
is promising, despite the critical current disorder inherent in their
fabrication process. Our study highlights the importance of theoretical
modeling in understanding the performance of superconducting devices and
provides insights that could inform the design of future SQIF arrays.

Suggested keywords

##  I Introduction

When fabricating 2D superconducting quantum interference filter (SQIF) arrays
(or superconducting quantum interference (SQUID) arrays) from YBCO thin films,
the critical currents $I_{ck}$ and electric resistances $R_{k}$ of the
Josephson junctions (JJs) exhibit significant variability and thus have wide
distributions. This $I_{c}$-$R$-disorder in the arrays is thought to result in
a decrease in the magnetic-field-to-voltage transfer function of the array,
which is an undesirable outcome as the goal is typically to achieve a transfer
function that is as large as possible for optimal performance in applications.
How detrimental such an $I_{c}$-$R$-disorder is in real devices, is still
unknown. This is in contrast to low-temperature superconducting thin-film
arrays, where the $I_{ck}$ and $R_{k}$ values can be well controlled.

In this paper, we aim to provide a quantitative analysis of the impact of the
$I_{c}$-$R$-disorder on the performance of two specific YBCO 2D SQIF arrays
that we have recently fabricated and measured. Our discussion will focus on
the effects of $I_{c}$-$R$-disorder in these particular arrays, while the
broader impact of $I_{c}$-$R$-disorder on 2D SQIF (or SQUID) arrays of varying
dimensions and with different screening and thermal noise parameters will be
addressed in a forthcoming publication.

Previously, we extensively modeled a structurally wide, thin-film 1D parallel
SQUID array with center current biasing and used a brute-force approach to
calculate the voltage-to-magnetic-field characteristics, as described in
reference [1]. In the present paper, we use an improved version of our
previous model, where we separate the superconducting current density into its
Meissner current, circulating loop currents, and bias-current components. This
separation allows us to save significant computational time and accurately
calculate larger 2D arrays with strong $I_{c}$-$R$-disorder.

Previous theoretical studies on 1D and 2D SQUID [2, 3, 4, 5], SQIF [6] and JJ
[7] arrays have often not considered the effects of thermal noise. Recent
research has shown that properly accounting for thermal noise is essential for
accurately calculating the voltage-to-magnetic-field characteristics and
transfer function of these arrays [8, 9, 10]. Only a few studies so far have
investigated the consequences of $I_{c}$-$R$-disorder in arrays [11, 5, 12],
but direct comparisons with experimental data from genuine arrays with wide
thin-film structures are lacking.

In this paper, we begin by discussing the layout, fabrication, and measurement
setup of our 2D SQIF arrays in Section II. We then present our experimental
results in Section III. In Section IV, we outline our theoretical model,
including the JJ phase dynamics and the stream function approach. In Section
V, we compare our experimental data with our model calculations and discuss
the insights gained from our detailed theoretical calculations with emphasis
on the $I_{c}$-$R$-disorder. We summarize our findings in Section VI.
Appendices are included to provide guidance on how to obtain the effective
hole areas and the geometric and kinetic inductances.

##  II 2D SQIF array layout, fabrication and measurement setup

Two different SQIF arrays were fabricated lithographically from a thin-film of
YBCO grown by e-beam evaporation on a 1 cm2 MgO substrate. The thickness of
the film was $d$ =113 nm and its superconducting transition temperature was
$T_{c}=85.9$ K. Before film deposition, steps were etched into the MgO surface
(see Fig. 1) using a well established technique based on argon milling [13,
14]. During YBCO thin-film deposition, grain boundaries form along the edges
of the sharp MgO steps, producing JJs. Films were then lithographically
patterned into two different 2D SQIF arrays referred to in the following as
arrays S1 and S2, which are displayed in Fig. 1.

Figure 1: 2D SQIF arrays S1 and S2. Etched rectangular regions in the
substrate are visible. Each of these regions has a pronounced sharp step edge
at the top to promote grain boundary thin-film growth and a very gradual step
on the bottom to avoid grain boundary formation in the busbar region.

Each array consists of $N_{s}=167$ identical 1D parallel SQIF arrays (the
rows) which are connected in series via superconducting YBCO thin-film
busbars. Each 1D SQIF array has $N_{p}$ = 6 JJs in parallel. Table I lists the
design parameters of S1 and S2, i.e. the SQIF hole height, the hole widths
(left to right), the busbar height, and the number $N_{s}$ of identical 1D
parallel SQIF arrays connected in series. The tracks leading to the JJ’s have
the same width as the JJs themselves (Table I). Both arrays have geometric
left-right symmetry. Because of the difficulty controlling grain boundary
growth, the values of the critical currents $I_{ck}$ and resistances $R_{k}$
($k$ = 1 to $N_{s}N_{p}$) of the JJs can vary widely [15]. This disorder in
$I_{ck}$ and $R_{k}$ negatively affects the applied magnetic field to voltage
transduction. In a 2D array, it is not possible to measure the $I_{ck}$ and
$R_{k}$ values individually and only consequences of the $I_{c}$-$R$-disorder
can be detected.

Table 1: SQIF array design parameters # | hole height | hole widths L-R | busbar height | JJ width | $N_{s}$
---|---|---|---|---|---
| [ $\mu$m ] | [ $\mu$m ] | [ $\mu$m ] | [$\mu$m] |
S1 | 15 | 1, 5, 13, 5, 1 | 16 | 2 | 167
S2 | 15 | 8, 12, 20, 12, 8 | 16 | 2 | 167

For measurements, the arrays were placed on a measurement probe and dipped
into a dewar of liquid nitrogen and zero field cooled from room temperature
down to 77 K. To screen out the earth’s magnetic field, the dewar was
surrounded by five layers of mu-metal shielding. The standard four-terminal
method was used to measure the DC voltages $V$ appearing between the ends
(between top and bottom) of the 2D SQIF arrays, while a bias current
$I_{b}^{tot}$ was injected at the top and exited the array at the bottom. A
solenoid surrounding the probe enabled us to generate a homogeneous
perpendicular applied magnetic field $B_{a}$ at the array. Both $V(B_{a})$ and
$V(I_{b}^{tot})$ of S1 and S2 were measured at 77 K.

##  III Experimental results

Figure 2 shows the experimental results of voltage $V$ versus the applied
magnetic field $B_{a}$ for arrays S1 and S2, where $B_{a}$ was varied from
$-50\,\mu\text{T}$ to $+50\,\mu\text{T}$.

Figure 2: Experimental data $V(B_{a})$ of arrays S1 and S2 for optimal bias
current $I_{b,opt}^{tot}$ (see Table II).

The bias currents $I_{b}^{tot}$ for S1 and S2 were chosen to optimize the
voltage modulation depths $\Delta V$. $\Delta V$ can be used as a measure of
the maximum transfer function $dV/dB_{a}|_{max}$ of the center dip [8]. The
values of the optimal bias current $I_{b,opt}^{tot}$ and $\Delta V$ are given
in Table II. Both $V(B_{a})$ curves are nearly reflection symmetric with
respect to the center dip. Because the rectangular holes in S2 are overall
larger than those in S1, the $V(B_{a})$ characteristics of array S2 shows more
oscillations in Fig. 2 than S1. To understand the origin of the complex shapes
of these two $V(B_{a})$ characteristics, we have performed theoretical
modeling which is outlined in detail below.

Table 2: SQIF array electrical parameters # |  $I_{b,opt}^{tot}$ [ $\mu$A ] |  $\Delta V$ [mV] |  $R^{tot}$ [ $\Omega$ ]
---|---|---|---
S1 | 76.474 | 8.52 | 406.0
S2 | 78.215 | 5.76 | 399.1

Figure 3 displays the results of our $V(I_{b}^{tot})$ characteristics
measurements in zero magnetic field for S1 and S2. The differential resistance
at large $I_{b}^{tot}$ defines the total resistance $R^{tot}$ of the array
which is given in Table II. Our theoretical modeling results in Section V will
give an explanation for the shape of $V(I_{b}^{tot})$.

Figure 3: Experimental data of voltage $V$ versus bias current $I_{b}^{tot}$
of S1 and S2 in zero applied magnetic field.

Section IV below will outline the theoretical model that we use to describe 2D
SQIF (and SQUID) arrays. We will see how our model allows us to extract from
our experimental data information about the JJ $I_{c}$-$R$-disorder in our
fabricated arrays. Furthermore, our model predicts how the voltage modulation
depth $\Delta V$, which can be used as a measure for the maximum transfer
function $dV/dB_{a}|_{max}$, decreases with increasing London penetration
depth $\lambda$ and increasing JJ $I_{c}$-$R$-disorder. Our model also
elucidates the role of the row-row mutual inductive coupling in 2D SQIF arrays
with busbars.

##  IV Theoretical model

###  IV.1 JJ phase dynamics

Our fabricated step-edge JJs have negligibly small capacitances and behave
like short JJs as the JJ width is much less than their Josephson penetration
depths at 77 K. Therefore, we use the resistively shunted junction model (RSJ
model) to describe our 2D SQIF arrays. The current $I_{k}$ flowing through the
$k^{th}$ JJ (Fig. 4), where $k$ = 1 to $N_{p}N_{s}$ (numbered from left to
right through the rows from top to bottom), is the sum of three different
currents: (i) the current through the intrinsic shunt resistance $R_{k}$, (ii)
the Josephson current $I_{ck}\sin\varphi_{k}$, and (iii) the noise current
$I_{Nk}$. Here $I_{ck}$ is the JJ critical current and $\varphi_{k}$ is the
gauge invariant phase difference. The noise current $I_{Nk}$ is due to the
Johnson noise voltage which appears across the JJ intrinsic shunt resistors
$R_{k}$ at finite temperature $T$.

Figure 4: Current flow schematic for a 2D SQIF array with uniform bias current
injection. The crosses mark the resistively shunted Josephson junctions. Here,
$N_{p}=6$ and $N_{s}=4$. For simplicity, only one circulating current $J_{m}$
and one current $I_{k}$ flowing through the $k^{th}$ JJ are indicated.

Employing the Josephson equation for the voltage across the $k^{th}$ JJ, one
obtains

$\frac{\Phi_{0}}{2\pi
R_{k}}\frac{d\varphi_{k}(t)}{dt}+I_{ck}\sin\varphi_{k}(t)+I_{Nk}(t)=I_{k}(t)\;.$
(1)

Here $\Phi_{0}$ is the flux quantum and $t$ the time. Introducing the
dimensionless time $\tau=2\pi RI_{c}t/\Phi_{0}$ where $R$ is the average JJ
resistance and $I_{c}$ the average critical JJ current, one can write Eq. (1)
in vector notation as

$\hat{\xi}^{-1}\frac{d\vec{\varphi}(\tau)}{d\tau}+\hat{\eta}\;\overrightarrow{\sin\varphi(\tau)}+\vec{i}_{N}(\tau)=\vec{i}(\tau)\;.$
(2)

Here $\hat{\xi}$ and $\hat{\eta}$ are $N_{s}Np\times N_{s}N_{p}$ diagonal
matrices with diagonal elements $\hat{\xi}_{kk}=R_{k}/R$ and
$\hat{\eta}_{kk}=I_{ck}/I_{c}$, respectively. The components of
$\vec{\varphi}$ are $\varphi_{k}$ and those of $\overrightarrow{\sin\varphi}$
are $\sin\varphi_{k}$, while the components of $\vec{i}_{N}$ are
$I_{Nk}/I_{c}$ and those of $\vec{i}$ are $I_{k}/I_{c}$.

From Kirchhoff’s law one obtains for a 2D array with uniform bias current
injection as indicated schematically in Fig. 4

$\vec{i}(\tau)=\hat{K}\,\vec{J}(\tau)/I_{c}+\vec{i}_{b}\;.$ (3)

Here $\vec{J}$ is the circulating current vector with components $J_{m}$
flowing around hole $m$ (Fig. 4), where $m$ = 1 to $N_{s}(N_{p}-1)$, numbered
from left to right through the rows from top to bottom. The transport bias
current vector $\vec{i}_{b}$ is $N_{s}N_{p}$ dimensional with identical
components $I_{b}/I_{c}$. The total bias current is $I_{b}^{tot}=N_{p}\,I_{b}$
where the time independent bias current is uniformly injected from the top as
shown in Fig 4. Because in our arrays $N_{s}>>N_{p}$, any effects arising from
the top and bottom bias current leads are neglected. The matrix $\hat{K}$ in
Eq. (3) is an $N_{s}N_{p}\times N_{s}(N_{p}-1)$ matrix of the form

$\hat{K}=\hat{I}_{N_{s}\times N_{s}}\otimes\hat{\kappa}\;,$ (4)

where $\otimes$ is a tensor product with $\hat{I}_{N_{s}\times N_{s}}$ an
$N_{s}\times N_{s}$ identity matrix and $\hat{\kappa}$ an
$N_{p}\times(N_{p}-1)$ matrix with elements
$(\hat{\kappa})_{ij}=\delta_{ij}-\delta_{i-1,j}$. Here $\delta_{ij}$ is the
Kronecker $\delta$.

Using the second Ginzburg-Landau equation one can relate neighboring JJ phases
within each row via the expression

$\frac{\Phi_{0}}{2\pi}(\varphi_{k+1}-\varphi_{k})=\Phi_{m}+\mu_{0}\,\lambda^{2}\oint_{C_{m}}\vec{j}\,d\vec{l}\;.$
(5)

Here $\lambda$ is the London penetration depth. The index $m$ is the index for
the circulating currents $J_{m}$ and for the holes. (The index $m$ is related
to the JJ index $k$ via the relation
$k=m+\text{Int}\\{m/(N_{p}-1)-\epsilon\\}$ where $\text{Int}\\{\\}$ is the
integer part of its argument and $\epsilon$ is an infinitesimally small
positive real number.) $\Phi_{m}$ in Eq. (5) is the magnetic flux through hole
$m$, and $\vec{j}$ is the super-current density in the array. The line
integration in Eq. (5) is performed counterclockwise along the contour $C_{m}$
which runs along the edges of the hole, and $\mu_{0}$ is the permeability of
free space.

The magnetic flux $\Phi_{m}$ in Eq. (5) is made up of four contributions,

$\Phi_{m}=\,\Phi^{(\text{app})}_{m}+\Phi^{(\text{Mei})}_{m}+\Phi^{(J)}_{m}+\Phi^{(I_{b})}_{m}\;.$
(6)

Here, $\Phi^{(\text{app)}}_{m}$ is the externally applied flux through hole
$m$, with $\Phi^{(\text{app})}_{m}=A_{m}B_{a}$, where $A_{m}$ is the hole area
and $B_{a}$ the applied homogeneous perpendicular magnetic field.
$\Phi_{m}^{\text{(Mei)}}$ is the flux through hole $m$ created by all the
Meissner shielding currents induced in response to the applied field $B_{a}$.
$\Phi_{m}^{(J)}$ is the flux produced in hole $m$ by all the circulating
currents $\vec{J}$, and $\Phi_{m}^{(I_{b})}$ is the flux through hole $m$
produced by the transport bias currents flowing vertically through all the
JJ’s from the top to the bottom (Eq. (3) and Fig. 4).

In a similar way, the supercurrent density $\vec{j}$ in Eq. (5) is made up of
three parts,

$\vec{j}\,=\,\vec{j}^{(Mei)}+\vec{j}^{(J)}+\vec{j}^{(I_{b})}\;.$ (7)

Here $\vec{j}^{(\text{Mei})}$ is the Meissner shielding current density,
$\vec{j}^{(J)}$ the current density from all the circulating currents, and
$\vec{j}^{(I_{b})}$ the transport bias current density in the array. The bias
current $I_{b}$ not only flows through the leads (Fig. 4), but also through
each JJ according to Eq. (3).

By using Eqs. (6) and (7), one can rewrite Eq. (5) as

$\frac{\Phi_{0}}{2\pi}(\varphi_{k+1}-\varphi_{k})=\Phi_{1,m}+\Phi_{2,m}+\Phi_{3,m}\;,$
(8)

where the fluxoids $\Phi_{1,m}$, $\Phi_{2,m}$ and $\Phi_{3,m}$ are defined as

$\Phi_{1,m}=\Phi_{m}^{(\text{app})}+\Phi_{m}^{(\text{Mei})}+\mu_{0}\lambda^{2}\oint_{C_{m}}\vec{j}^{(\text{Mei})}\,d\vec{l}\;,$
(9)
$\Phi_{2,m}=\Phi_{m}^{(\text{J})}+\mu_{0}\lambda^{2}\oint_{C_{m}}\vec{j}^{(\text{J})}\,d\vec{l}$
(10)

and

$\Phi_{3,m}=\Phi_{m}^{(I_{b})}+\mu_{0}\lambda^{2}\oint_{C_{m}}\vec{j}^{(I_{b})}\,d\vec{l}\;.$
(11)

The fluxoid $\Phi_{1,m}$ in Eq. (9) can be written as

$\Phi_{1,m}=B_{a}A_{m}^{\text{eff}}\;,$ (12)

where $A_{m}^{\text{eff}}$ is the effective hole area of hole $m$, with

$A_{m}^{\text{eff}}=\frac{B_{a}A_{m}+\Phi^{(\text{Mei})}+\mu_{0}\lambda^{2}\oint_{C_{m}}\vec{j}^{(\text{Mei})}\,d\vec{l}}{B_{a}}\,.$
(13)

The effective hole area $A_{m}^{\text{eff}}$ does not depend on $B_{a}$,
because both $\Phi^{(\text{Mei})}$ and $\vec{j}^{(\text{Mei})}$ are
proportional to $B_{a}$. How we determine $A_{m}^{\text{eff}}$ is outlined
further below.

Furthermore, the fluxoid $\Phi_{2,m}$ in Eq. (10) can be written as

$\Phi_{2,m}=(\hat{L}\,\vec{J})_{m}\;,$ (14)

where $\hat{L}$ is the $N_{s}(N_{p}-1)\times N_{s}(N_{p}-1)$ inductance square
matrix of the array, with each matrix element the sum of a geometric and
kinetic term. How we determine the inductance matrix $\hat{L}$ is outlined
further below.

And, the fluxoid $\Phi_{3,m}$ in Eq. (11) can be written as

$\Phi_{3,m}=(\vec{L}^{I_{b}})_{m}\,I_{b}\;,$ (15)

where we name $\vec{L}^{I_{b}}$ the bias current inductance vector which has
$N_{s}(N_{p}-1)$ components. Each component is the sum of a geometric and
kinetic term. How we determine the bias current inductance vector
$\vec{L}^{I_{b}}$ is outlined further below.

Finally, using Eqs. (2), (3), and (8) - (15), we derive the following system
of nonlinear stochastic differential equation of first order for the time
evolution of the JJ phase difference vector $\vec{\varphi}(\tau)$

$\hat{\xi}^{-1}\frac{d\vec{\varphi}(\tau)}{d\tau}+\hat{\eta}\;\overrightarrow{\sin\varphi(\tau)}+\vec{i}_{N}(\tau)=$
$\vec{i_{b}}\,+\,\frac{\hat{K}\hat{L}^{-1}}{I_{c}}\,[\;\frac{\Phi_{0}}{2\pi}\hat{D}\vec{\varphi}(\tau)-B_{a}\vec{A}^{\text{eff}}-\vec{L}^{I_{b}}\,I_{b}\;]\;.$
(16)

Here $\hat{D}$ is an $N_{s}N_{p}\times N_{s}(N_{p}-1)$ matrix of the form

$\hat{D}=\hat{I}_{N_{s}\times N_{s}}\otimes\hat{\delta}\;,$ (17)

where $\hat{\delta}$ is an $(N_{p}-1)\times N_{p}$ matrix with elements
$(\hat{\delta})_{ij}=-\delta_{ij}+\delta_{i,j-1}$. The differential equation
Eq. (16) is the key equation of this paper.

The set of dynamic equations for the JJ phase differences $\varphi_{k}$ (Eq.
(16)) can also be written with the help of a so-called washboard potential
$U(\\{\varphi_{k}\\})$ as [7]

$\frac{\Phi_{0}^{2}}{4\pi^{2}R_{k}}\frac{d\varphi_{k}}{dt}=-\frac{\partial
U(\\{\varphi_{k}\\})}{\partial\varphi_{k}}\;,$ (18)

where

$U(\\{\varphi_{k}\\})=-\sum_{k=1}^{N_{s}N_{p}}[\,\frac{\Phi_{0}}{2\pi}I_{ck}\cos\varphi_{k}+\frac{\Phi_{0}}{2\pi}(I_{b}-I_{Nk})\varphi_{k}\,]$
(19)
$+\,\frac{1}{2}\sum_{m=1,m^{\prime}=1}^{N_{s}(N_{p}-1)}M_{m}\,\hat{L}^{-1}_{m,m^{\prime}}\,M_{m^{\prime}}\;,$

with

$M_{m}=\frac{\Phi_{0}}{2\pi}(\varphi_{k+1}-\varphi_{k})-A_{m}^{\text{eff}}B_{a}-L_{m}^{I_{b}}I_{b}\;,$
(20)

where $k$ is a function of $m$ as stated above.

The potential $U$ can be obtained from energy considerations or conjectured
directly from Eq. (16). Evaluating the right hand side of Eq. (18) leads to
Eq. (16).

The time averaged voltage $V$ of the SQIF array, measured between its top and
bottom ends, is given by

$V=R\,I_{c}\lim_{{\tau}\rightarrow\infty}\,\frac{1}{\tau}\,\frac{1}{N_{p}}\sum_{k=1}^{N_{s}N_{p}}[\,\varphi_{k}(\tau)-\varphi_{k}(0)\,]\;.$
(21)

The noise current vector $\vec{i}_{N}(\tau)$, which appears in Eq. (16), was
implemented using uncorrelated Gaussian random variables with the mean-square
deviations $2\,\Gamma_{k}/\Delta\tau$. Here $\Delta\tau$ is the time-step
interval used in the finite difference method used to numerically solves Eq.
(16), and $\Gamma_{k}$ is the noise strength defined as

$\Gamma_{k}=\frac{R}{R_{k}}\,\frac{2\pi k_{B}T}{I_{c}\Phi_{0}}\,,$ (22)

where $k_{B}$ is the Boltzmann constant and $T$ the operating temperature of
the array.

###  IV.2 Stream function approach

In this section we show how, for superconducting thin-film arrays with
extended busbar structures (Fig. 1), one can use the stream function approach
[16, 17, 18, 19, 1, 20] to calculate the effective hole areas
$\vec{A}^{\text{eff}}$ (Eq. (13)), the array inductance matrix $\hat{L}$ (Eq.
(14)), and the bias current inductance vector $\vec{L}^{I_{b}}$ (Eq. (15)).
These three physical quantities are needed in order to numerically solve the
set of differential equations in Eq. (16) and to obtain the time evolution of
the JJ phase differences $\varphi_{k}(\tau)$.

For a superconducting thin-film of thickness $d$ in the $xy$-plane, the super-
current density $\vec{j}(x,y,z)$ becomes independent of the perpendicular
direction $z$ if $\lambda>d$ [21, 18]. In our case we have $d=0.113\,\mu$m and
$\lambda=0.42\,\mu$m (as shown later) and thus $\lambda>d$. In this case one
can introduce a 2D stream function $g(x,y)$ which is defined via the $x$ and
$y$ components of $\vec{j}$ as

$j_{x}=\frac{1}{d}\frac{\partial g}{\partial
y}\;\;\text{and}\;\;j_{y}=-\frac{1}{d}\frac{\partial g}{\partial x}\;.$ (23)

Using Eq. (23) together with the second London equation and Biot-Savart’s law
in 2D, one obtains by applying partial integration a second-order linear
Fredholm intregro-differential equation for the stream function $g(x,y)$ of
the form

$\frac{\lambda^{2}}{d}\,(\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial y^{2}})\,g(x,y)\;+$
$\frac{1}{4\pi}\int_{\Omega}Q(x,y,x^{\prime},y^{\prime})\;g(x^{\prime},y^{\prime})\;dx^{\prime}dy^{\prime}-f_{s}(x,y)\,=\,\frac{B_{a}}{\mu_{0}}\;,$
(24)

where $B_{a}$ is the applied magnetic field pointing in $z$ direction and

$f_{s}(x,y)\;=$ (25)
$\frac{1}{4\pi}\oint_{\partial\Omega}\frac{g(x^{\prime},y^{\prime})}{\sqrt{(x-x^{\prime})^{2}+(y-y^{\prime})^{2}}^{\;3}}\,\left({\begin{array}[]{cc}x-x^{\prime}\\\
y-y^{\prime}\\\ \end{array}}\right)\vec{n}\;dl^{\prime}\,.$

In Eq. (24), $\Omega$ is the integration domain, i.e. the superconducting thin
film area, and $\partial\Omega$ is the domain boundary which includes the
edges of the holes. The vector $\vec{n}$ in Eq. (25) is a normal vector in the
$xy$-plane, perpendicular to the boundary, and $\vec{n}$ points outwards, away
from the domain $\Omega$. In Eq. (25), $dl^{\prime}$ is a line element.

The kernel $Q$ in Eq. (24) can be written as [1]

$Q(x,y,x^{\prime},y^{\prime})=$ (26) $\frac{1}{\Delta x\Delta
y}\,\Bigg{[}\,\Big{[}\frac{\sqrt{\bar{x}^{2}+\bar{y}^{2}}}{\bar{x}\,\bar{y}}\Big{]}_{x^{\prime}-x-\Delta{x}/2}^{x^{\prime}-x+\Delta{x}/2}(\bar{x})\Bigg{]}_{y^{\prime}-y-\Delta{y}/2}^{y^{\prime}-y+\Delta{y}/2}(\bar{y})\;\;,$

with the definition

$\Bigg{[}\,\Big{[}f(x,y)\Big{]}_{s_{1}}^{s_{2}}(x)\Bigg{]}_{s_{3}}^{s_{4}}(y)\,=$
(27) $f(s_{2},s_{4})-f(s_{1},s_{4})-f(s_{2},s_{3})+f(s_{1},s_{3})\;,$

where $\Delta x$ and $\Delta y$ are the grid spacings used when solving
solving Eq. (24) numerically. We used $\Delta x$ = $\Delta y$.

Having determined the stream function $g(x,y)$ allows us to calculate the
magnetic fluxes $\Phi_{m}$ of the holes $m$ as well as the
$\oint_{\partial\Omega_{m}}\vec{j}\,d\vec{l}$ terms. One obtains

$\Phi_{m}=\mu_{0}\int_{\Omega_{m}}h(x,y)\;dx\,dy\;+\,\Phi_{m}^{(\text{app})},$
(28)

where $h(x,y)$ is defined as

$h(x,y)=f_{s}(x,y)\,-\,\frac{1}{4\pi}\int_{\Omega}\,Q(x,y,x^{\prime},y^{\prime})\;g(x^{\prime},y^{\prime})\;dx^{\prime}dy^{\prime}\;,$
(29)

and

$\oint_{\partial\Omega_{m}}\vec{j}\,d\vec{l}\,=\,\frac{1}{d}\,\oint_{\partial\Omega_{m}}\,(\,\frac{\partial
g(x,y)}{\partial y}\,dx\,-\,\frac{\partial g(x,y)}{\partial x}\,dy\,)\;.$ (30)

By choosing different boundary conditions for the stream function $g(x,y)$,
one can separately address the three different current densities in Eq. (7),
i.e. $j^{(Mei)}$, ${j^{(J)}}$ and ${j^{(I_{b})}}$, and then obtain
$\vec{A}^{\text{eff}}$, $\hat{L}$ and $\vec{L}^{I_{b}}$ from Eqs. (13) - (15).
Further details about how to obtain $\vec{A}^{\text{eff}}$, $\hat{L}$ and
$\vec{L}^{I_{b}}$ are given in Appendices A, B and C.

##  V Model versus experimental results

There are $2N_{s}N_{p}+1$ unknown parameters in our model. These unknown
parameters are the $N_{s}N_{p}=167\times 6=1002$ values for the JJ critical
currents $I_{ck}$, the $N_{s}N_{p}=1002$ values for the JJ resistances
$R_{k}$, and the London penetration depth $\lambda$. Experimentally it has
been found that in YBCO JJs, $I_{ck}$ and $R_{k}$ are approximately anti-
correlated according to the empirical law $I_{ck}R_{k}\sim j_{ck}^{1/2}$,
where $j_{ck}$ is the JJ critical current density [22, 14]. By exploiting this
empirical law, we obtain for the diagonal matrix $\hat{\xi}$ in Eq. (16) the
matrix elements $\hat{\xi}_{kk}=\hat{\eta}_{kk}^{-1/2}$. This reduces the
number of unknown model parameters by almost a factor of 2, from
$2\,N_{s}N_{p}+1$ to $N_{s}N_{p}+1$.

The critical currents of YBCO step-edge junctions have a spread in values and
previous measurements have found wide $I_{c}$-distributions of approximately
Gaussian or even log-normal [15] shape. Therefore, when solving Eq. (16), we
treat the $N_{s}N_{p}$ matrix elements $\hat{\eta}_{kk}=I_{ck}/I_{c}$ as log-
normal random variables. A log-normal distribution $p(\eta)$ with mean value 1
and standard deviation $\sigma$ has the form

$p(\eta)=\frac{1}{\eta\,\gamma\,\sqrt{2\pi}}\,\exp\frac{-(\ln\eta-\mu)^{2}}{2\,\gamma^{2}}\;.$
(31)

Here $\mu=-\,\gamma^{2}/2$ and $\gamma=\sqrt{\ln(1+\sigma^{2})}$. A log-normal
distribution is similar to a Gaussian distribution for small $\sigma$. In our
modeling, we cannot use Gaussian random variables for $\hat{\eta}_{kk}$
because they would lead to unphysical negative $I_{ck}/Ic$ values,
particularly for large $\sigma$. In the following we call the standard
deviation $\sigma$ the $I_{c}$-spread.

The number of unknown parameters is reduced slightly further by the relation
between the measured total resistance $R^{tot}$ (Table II) and the unknown
average JJ resistance $R$. This relation is given by

$R=\frac{R^{tot}}{\sum_{i=1}^{N_{s}}\,(\sum_{j=1}^{Np}\hat{\eta}_{kk\leftarrow
ij}^{1/2})^{-1}}\;\,.$ (32)

Here the symbol $kk\leftarrow ij$ means that the $ij$ indices (where $i=1$ to
$N_{s}$ is the row index and $j=1$ to $N_{p}$, the column index) are mapped
onto the $k$ index which runs from 1 to $N_{s}N_{p}$. For arrays with
sufficiently large $N_{s}N_{p}$, due to the introduction of $p(\eta)$ and Eq.
(32), one no longer is dealing with $N_{s}N_{p}+1$ parameters but instead one
only has to deal with three parameters, which are $I_{c}$, $\lambda$ and
$\sigma$. Because S1 and S2 were produced from the same YBCO thin film on the
same substrate, one can expect $\lambda$ and $\sigma$ to be the same for S1
and S2. But because of the finite values $N_{p}=6$ and $N_{s}=167$, one can
expect a slight difference between the average $I_{c}$ values of array S1 and
S2. The values of $I_{c}$, $\lambda$ and $\sigma$ for S1 and S2, that were
found to best model our experimental data, are listed in Table III. These
values were found by probing many different $I_{c}$, $\lambda$ and $\sigma$
combinations and comparing by eye the calculated $V(B_{a})$ and
$V(I_{b}^{tot})$ curves with the experimental data curves. The search for the
$I_{c}$, $\lambda$ and $\sigma$ parameters was made easier by exploiting the
facts that the $V(I_{b}^{tot})$ curves are very sensitive to $I_{c}$ and that
decreasing $\lambda$ shifts $V(B_{a})$ upwards and increases the voltage
modulation depth $\Delta V$, whereas increasing $\sigma$ only decreases
$\Delta V$. Because of the $I_{c}$-disorder and the finite size of the arrays,
the parameters $I_{c}$, $\lambda$ and $\sigma$ can only be determined with our
procedure with some uncertainties which we estimated to be about 2%, 5% and
10%, respectively. Previously, in DC-SQUIDs [23] that were made from similar
YBCO thin films, a value of $\lambda(77\text{K})\simeq 0.39\,\mu$m was found,
which is close to $\lambda=0.42\,\mu$m used here.

Table 3: Best fit parameters # |  $I_{c}$ [$\mu$A] |  $\lambda$ [$\mu$m] |  $\sigma$ [$\mu$m]
---|---|---|---
S1 | 19.3 | 0.42 | 0.5
S2 | 20.2 | 0.42 | 0.5

While doing the calculations we found that the inductance matrix elements
between rows are quite small because the superconducting busbars sufficiently
separate the holes of adjacent 1D parallel SQIF arrays (rows). We noticed that
neglecting the small row-row mutual inductances made no noticeable difference
to any of our results. Therefore, our 2D arrays behave like $N_{s}$
independent 1D parallel SQIF arrays where the 2D array DC voltage is simply
the sum over all the row DC voltages. Being allowed to neglect the row-row
coupling simplifies Eq. (16) and dramatically reduced the computational time
needed to calculate the time evolution of the JJ phases.

Though each of our 2D SQIF arrays contains $N_{s}N_{p}=1002$ JJs, $N_{s}$ is
still not large enough for $V(B_{a})$ and $V(I_{b}^{tot})$ to become
independent of the chosen $I_{c}$-disorder set {$\hat{\eta}_{kk}$}. We have
investigated 40 different randomly chosen sets {$\hat{\eta}_{kk}$} for both S1
and S2, and selected the two $I_{c}$-disorder sets that fitted our S1 and S2
experimental data best.

Figure 5 shows calculated $V(B_{a})$ characteristics for S1 and S2, using the
parameters in Table III, for 10 different randomly chosen $I_{c}$-disorder
sets {$\hat{\eta}_{kk}$} for each array. The $V(B_{a})$ curves vary, as each
$I_{c}$-disorder set {$\hat{\eta}_{kk}$} produces a somewhat different result.
Figure 6 displays the $I_{c}$-disorder set {$\hat{\eta}_{kk}$} used to
calculate the dotted $V(B_{a})$ characteristics of S1 in Fig. 5. Due to the
long tail of a log-normal distribution with $\sigma=0.5$, values of
$\hat{\eta}_{kk}$ as large as 4 can be seen.

Figure 5: $V(B_{a})$ characteristics of arrays S1 and S2 for 10 different
$I_{c}$-disorder sets {$\hat{\eta}_{kk}$}. Figure 6: $I_{c}$-disorder set
{$\hat{\eta}_{kk}=I_{ck}/I_{c}$} in array S1 used for calculating the dotted
curve in Fig. (5). The array has $N_{p}=6$ JJ columns and $N_{s}=167$ JJ rows.

The calculated dotted curves from Fig. 5 for S1 and S2, together with our
experimental S1 and S2 data from Fig. 2, are displayed in Fig. 7. Figure 7
shows a key result of our paper. The quantitative agreements between
experimental data and calculation for both S1 and S2 are exceptional. Still,
we could have improved the agreement even further by searching for more
optimal $I_{c}$-disorder sets {$\hat{\eta}_{kk}$} . It turns out that the
effective hole areas $A_{m}^{\text{eff}}$ are important ingredients to
describe the experimental $V(B_{a})$ data correctly. As can be seen in Fig. 7,
the calculated curves are slightly more stretched along the $B_{a}$ axis than
the experimental data. After the calculation, we found out that this
discrepancy can be partially explained by the fact that the layout size of the
arrays was actually about 2% larger than the values given in Table I.
Furthermore, our calculations showed that the term in Eq. (16) which contains
the bias current inductance $\vec{L}^{I_{b}}$ can be neglected without
affecting our calculated $V(B_{a})$ characteristics.

Figure 7: Comparison between experimental and calculated $V(B_{a})$
characteristics for array S1 and S2.

To gain a better understanding of the role of the London penetration depth
$\lambda$, Fig. 8 shows, for S1 and S2, how $V(B_{a})$ (for positive $B_{a}$)
changes with $\lambda$ while all other parameters stay unchanged. Increasing
$\lambda$ causes the voltage modulation depth $\Delta V=V_{max}-V_{min}$ to
decrease. This is mainly due to the increase in the kinetic inductances which
are proportional to $\lambda^{2}$. This is similar to a DC-SQUID where an
increase of the screening parameter $\beta_{L}$ leads to a decrease of $\Delta
V$ [24].

Figure 8: Calculated $V(B_{a})$ characteristics of array S1 and S2 for three
different values of $\lambda$.

Figure 9 shows how rapidly the voltage modulation depth $\Delta V(\lambda)$
for S1 and S2, which is a measure of the maximum transfer function
$dV/dB_{a}|_{max}$ [8], decreases with increasing $\lambda$. Note that
$\lambda$ enters our calculations via the factor $\lambda^{2}/d$ (Eq. (24),
and Eq. (30) with Eq. (5)), which is known as the Pearl penetration depth
[25]. Thus, increasing the film thickness $d$, as long as $\lambda>d$ [21,
18], would further increase the modulation depth $\Delta V$ of our arrays,
assuming that all other parameters stay unchanged. Figure 10 shows, for S1 and
S2, how the $V(B_{a})$ characteristics (for positive $B_{a}$) change with the
$I_{c}$-spread $\sigma$ while all other parameters stay unchanged. Increasing
$\sigma$ smoothes the $V(B_{a})$ curves, which decreases the voltage
modulation depth $\Delta V$. The underlying cause

Figure 9: Calculated voltage modulation depth $\Delta V$ versus $\lambda$ of
array S1 and S2.

Figure 10: Calculated $V(B_{a})$ characteristics of array S1 and S2 for three
different $I_{c}$-spreads $\sigma$.

is revealed by looking at the $V(\Phi^{\text{(app)}})$ characteristics of the
simple DC-SQUID ($N_{s}=1,N_{p}=2$) with $R_{1}\neq R_{2}$ and $I_{c1}\neq
I_{c2}$ [24, 26]. Here, the $R$-asymmetry causes a skewing of the
$V(\Phi^{\text{(app)}})$ response [24] and the $I_{c}$-asymmetry a
$\Delta\Phi^{\text{app}}$-shifting to the left or right by an amount
$\Delta\Phi^{\text{(app)}}=L_{s}\,|\,I_{c1}-I_{c2}\,|\,/\,2$ where $L_{s}$ is
the self-inductance of the DC-SQUID loop [26]. In addition to the reflection
symmetry breaking, the modulation depth decreases with increasing
$I_{c}$-asymmetry [24, 26]. The modulation depth is less affected by the
$R$-asymmetry. In our 2D SQIF arrays each row has 6 JJ’s in parallel and 5
holes. As the $I_{ck}$’s and $R_{k}$’s are all different, this produces in
each row a different reflection asymmetry, i.e. $V_{row}(B_{a})\neq
V_{row}(-B_{a})$, due to different skewing and $B_{a}$-shifting. Because of
the very weak inductive coupling between rows, as mentioned above, the time
averaged voltages $V_{row}(B_{a})$ across rows can simply be added up to give
the total array voltage $V(B_{a})$. Adding up a large number of different row
voltages $V_{row}(B_{a})$ causes a $V(B_{a})$ smoothing which reduces the
modulation depth $\Delta V$ of the array but also reduces the reflection
asymmetry.

Figure 11 displays $\Delta V(\sigma)$ normalized by $\Delta V(\sigma=0)$. The
figure shows that an $I_{c}$-spread $\sigma=0.5$ reduces $\Delta
V(\sigma)/\Delta V(0)$ in the case of S1 from 1 to about 0.55 and in the case
of S2 from 1 to about 0.4. Figure 11 reveals that $I_{c}$-$R$-disorder has a
stronger detrimental effect on $\Delta V$ in an array with larger inductances
like S2 than in an array with smaller inductances like S1. Figure 11 is one of
the key results of this paper.

Figure 11: The calculated normalized voltage modulation depth $\Delta
V(\sigma)\,/\,\Delta V(0)$ versus $I_{c}$-spread $\sigma$ for array S1 and S2.

In Appendix D we discuss the weak reflection asymmetry of $V(B_{a})$.

Figure 12 shows the experimental and calculated $V(I_{b}^{tot})$
characteristics in zero applied magnetic field for S1 an S2. The calculation
does not perfectly reproduce the experimentally observed $V(I_{b}^{tot})$ at
low $I_{b}^{tot}$. Better agreement can possibly be achieved by modifying the
log-normal distribution $p(\eta)$ in Eq. (31) such that low $I_{c}$ values
become more emphasized. The values of $N_{p}I_{c}$ for S1 and S2 are indicated
in Fig 12. In contrast to $V(B_{a})$, the calculated shape of $V(I_{b}^{tot})$
is not very sensitive to a particular choice of the $I_{c}$-disorder set
{$\hat{\eta}_{kk}$} for a given $\sigma$, but is very sensitive to the choice
of $I_{c}$, i.e. the average $I_{ck}$ value.

Figure 13 again displays $V(I_{b}^{tot})$ for S1 but here the calculated
additional dotted curves reveal how the parameters $T$ and $\sigma$ affect the
shape of $V(I_{b}^{tot})$. Going from curve (1) to (2) shows the effect of
temperature $T=77$ K, and going from (2) to (3) the effect of $\sigma=0.5$.

Figure 12: Experimental and calculated $V(I_{b}^{tot})$ characteristics of
array S1 and S2. Figure 13: Experimental and calculated $V(I_{b}^{tot})$
characteristics of array S1. The parameters used to calculate the 3 dotted
blue curves are (1) $T=\sigma=0$, (2) $T=77$ K, $\sigma=0$, and (3) $T=77$ K,
$\sigma=0.5$.

##  VI Summary and Conclusion

We measured the $V(B_{a})$ characteristics at optimal bias currents and the
$V(I_{b}^{tot})$ characteristics at $B_{a}=0$ of two 2D SQIF arrays with
$N_{s}=167$ and $N_{p}=6$.

To model our experimental data, we developed a theoretical model that takes
into account the Johnson noise in the JJ resistances $R_{k}$, the spreads in
$I_{ck}$ and $R_{k}$, and the wide, thin-film busbar structure of the arrays.
The model was made numerically efficient by dividing the superconducting
current density into its Meissner current, circulating current, and bias
current parts. Using the stream function approach for $\lambda>d$, we
calculated the effective areas of the array holes, inductance matrices and
bias current inductance vectors. The calculations revealed that the busbar
height was sufficiently large such that the inductive row-row coupling could
be neglected. This significantly simplified the numerical calculations and
made the investigation of the $I_{c}$-$R$-disorder of our $N_{s}=167$ arrays
computationally possible.

Empirical findings that $R_{k}$ can be expressed in terms of $I_{ck}$ reduced
the parameter space, and we could identify $\\{I_{ck}\\}$ sets from a log-
normal distribution that produced $V(B_{a})$ characteristics that agreed well
with our experimental data.

Our calculations demonstrated that the voltage modulation depth $\Delta V$,
which is a measure for the maximum transfer function, decreases with
increasing $\lambda$ due to the increase in kinetic inductances which are
proportional to $\lambda^{2}/d$.

Most of all, our calculations showed how $\Delta V$ decreases with increasing
$I_{c}$-spread $\sigma$, where the decrease was more rapid for the high
inductance array S2. Our experimental data were best described by
$\sigma=0.5$, showing that the $I_{c}$-$R$-disorder in our arrays reduced
$\Delta V$ by 45% in S1 and 60% in S2.

We observed slight reflection asymmetries in our experimental $V(B_{a})$
characteristics, caused by the $I_{c}$-$R$-disorder, and our model predicted a
similar degree of asymmetry as observed experimentally. Generally, the
refection asymmetry increases with increasing $I_{c}$-spread $\sigma$, but
decreases when $N_{s}$ is increased.

The calculated $V(I_{b}^{tot})$ characteristics did not fully agree with our
experimental data at low $I_{b}^{tot}$, possibly due to the choice of a log-
normal distribution that did not sufficiently emphasize low $I_{ck}$ values.
Modifying our log-normal distribution would give better agreement with
$V(I_{b}^{tot})$ without affecting much the $V(B_{a})$ results. Our
calculations also showed how the $V(I_{b}^{tot})$ characteristic changes with
temperature ($T=0$ to $T=77$ K) and with $\sigma$ ($\sigma=0$ to
$\sigma=0.5$).

Our findings suggest that the use of YBCO thin films in SQIF arrays is
promising, as the critical current disorder inherent in their fabrication
process reduces the transfer function in our case only by about a factor of 2.
This is one of our important results. Our study highlights the importance of
theoretical modeling in understanding the performance of superconducting
devices and provides insights that could inform the design of future SQIF
arrays. While our detailed theoretical modeling provided valuable insight into
the role of $I_{c}$-$R$-disorder for the specific arrays tested, more general
investigations are needed to fully understand its role as a function of
screening parameter $\beta_{L}$, noise strength $\Gamma$, kinetic inductance
fraction, $N_{s}$ and $N_{p}$ values, and different array planar geometries.

###### Acknowledgements.

The authors would like to thank Marc Gali Labarias for stimulating
discussions.

##  Appendix A Effective hole areas

By choosing any constant applied magnetic field, $B_{a}>0$, and the boundary
condition $g|_{\partial\Omega}=0$, for all boundaries $\partial\Omega$, one
can calculate numerically, by solving Eq. (24), the part of $g(x,y)$ that
corresponds to the Meissner shielding current density $j^{(Mei)}$. In this
case $f_{s}=0$ (Eq. (25)). Stream-lines of $j^{(Mei)}$ in array S1 and S2 are
shown in Fig. 14. Using Eqs. (28) and (30), the fluxoids $\Phi_{1,m}$ (Eq.
(9)) of the different holes are determined and the effective hole area
enhancement factors $A_{m}^{\text{eff}}/A_{m}$ (Eq. (13)) are calculated.

Figure 14: Meissner current streamlines in S1 and S2.

$A_{m}^{\text{eff}}/A_{m}$ for different holes $m$ are displayed in Fig. 15,
where the geometric and kinetic parts are indicated. Figure 15 shows that the
kinetic part plays an important role.

Figure 15: Effective hole area enhancement factor $A_{m}^{eff}/A_{m}$ of S1
and S2. Geometric and kinetic parts are indicated.

##  Appendix B Inductance matrix $\hat{L}$

One can determine the part of $g(x,y)$, which corresponds to a
counterclockwise circulating current $J_{m}$ around hole $m$, by choosing
$B_{a}=0$ and solving Eq. (24) numerically with the appropriate boundary
condition. The boundary condition is $g|_{\Omega_{m}}=J_{m}$ ($J_{m}>0$ but
arbitrary), and $g$ increases linearly from 0 to $J_{m}$ along the left JJ and
decreases linearly from $J_{m}$ to 0 along the right JJ. Everywhere else
$g|_{\partial\Omega}=0$. In this case, $f_{s}$ in Eq. (25) can be obtained
analytically as $f_{s}(x,y)=[\,P^{u}_{m}(x,y)+P^{u}_{m}(x,-y)\,]\,J_{m}$,
where the formulas for $P^{u}_{m}$ can be found in Appendix B of ref. [1].
Examples of current stream-lines of currents flowing around a hole $m=2$ in S1
and S2 are shown in Fig. 16. One can use Eqs. (28) and (30) to determine the
fluxoids $\Phi_{2,m}$ (Eq. (10)) generated in the different holes, and then
with Eq. (14) all the matrix elements of the inductance matrix $\hat{L}$ are
calculated. The matrix elements $\hat{L}_{ij}$ of S1 and S2, for the holes in
a single row, are plotted in Fig. 17 and their geometric and kinetic parts are
indicated. The large self-inductances are seen along the diagonal while the
negative mutual inductances are off-diagonal. The row-row mutual inductance
matrix elements are found to be very small and are not displayed here. Figure
17 shows that the kinetic part plays a dominant role.

Figure 16: Current streamlines in array S1 and S2 to determine the inductance
matrix $\hat{L}$.

Figure 17: Inductance matrix elements $\hat{L}_{ij}$ of a single row, $i=1$ to
5 and $j=1$ to 5, for S1 and S2. The geometric and kinetic inductance parts
are shown in different colours.

##  Appendix C Inductance vector $\vec{L}^{I_{b}}$

One can determine the part of $g(x,y)$ that corresponds to the bias transport
current density $\vec{j}^{(I_{b})}$ by choosing $B_{a}=0$ and solving Eq. (24)
numerically with the appropriate boundary condition. The boundary condition is
$g=\pm N_{p}I_{b}/2$ along the left/right outer edges of the array, and, along
the inner edges of the holes, $g=-N_{p}I_{b}/2+iI_{b}$ with $i=1$ to
$N_{p}-1$, while $g$ increases linearly along the JJ’s between neighboring
holes. Again, $f_{s}$ can be obtained analytically and can be expressed in
terms of $P^{u}_{m}(x,y)$ [1]. Streamlines of $\vec{j}^{(I_{b})}$ in S1 and S2
are shown in Fig. 18. Using Eqs. (28) and (30), the fluxoids $\Phi_{3,m}$ (Eq.
(11)) of the different holes can be determined, and with Eq. (15) the
components $L_{m}^{I_{b}}$ of the inductance vector $\vec{L}^{I_{b}}$ are
calculated.

Figure 18: Current streamlines of $\vec{j}^{(I_{b})}$ in S1 and S2 to
determine $\vec{L}^{I_{b}}$.

Figure 19: Components $L_{m}^{I_{b}}$ of the bias current inductance vector
$\vec{L}^{I_{b}}$. The geometric and kinetic parts are shown.

Figure 19 shows $L_{m}^{I_{b}}$ for the different holes $m$ along a row, and
the geometric and kinetic parts are indicated.

##  Appendix D Reflection asymmetry of $V(B_{a})$

Figure 20: Experimental $V(B_{a})$ and $V(-B_{a})$ of array S1 and S2 to show
the reflection asymmetry.

Figure 21: Calculated $V(B_{a})$ and $V(-B_{a})$ of array S1 and S2 to show
the reflection asymmetry.

Figure 20 shows the experimental reflection asymmetry of $V(B_{a})$ by
plotting $V(B_{a})$ and $V(-B_{a})$ for S1 and S2. Small differences between
$V(B_{a})$ and the reflected curve $V(-B_{a})$ are noticeable. Figure 21
displays the calculated $V(B_{a})$ and $V(-B_{a})$ curves. One can define the
degree $\delta$ of the reflection asymmetry as

$\delta=\frac{\int_{-\Delta B/2}^{\Delta
B/2}\,\left|\;V(B_{a})-V(-B_{a})\;\right|\,dB_{a}}{\Delta V\Delta B}\;,$ (33)

with the integration interval $\Delta B=100\,\mu$T. The experimental $\delta$
for S1 and S2 are $\delta=0.0137$ and 0.015, respectively, while the
calculated ones are $\delta=0.0127$ and 0.012. We found from our calculations
that $\delta$ increases when the $I_{c}$-spread $\sigma$ increases.

## References

  * Müller and Mitchell [2021] K.-H. Müller and E. E. Mitchell, Theoretical model for parallel squid arrays with fluxoid focusing, Phys. Rev. B 103, 054509 (2021).
  * Cybart  _et al._ [2012] S. A. Cybart, T. N. Dalichaouch, S. M. Wu, S. M. Anton, J. A. Drisko, J. M. Parker, B. D. Harteneck, and R. C. Dynes, Comparison of measurements and simulations of series-parallel incommensurate area superconducting quantum interference device arrays fabricated from YBa2Cu3O7-δ ion damage Josephson junctions, J. Appl. Phys. 112, 063911 (2012).
  * Dalichaouch  _et al._ [2014] T. N. Dalichaouch, S. A. Cybart, and R. C. Dynes, The effects of mutual inductances in two-dimensional arrays of Josephson junctions, Supercond. Sci. Technol. 27, 065006 (2014).
  * Taylor  _et al._ [2016] B. J. Taylor, S. A. E. Berggren, M. C. O’Brian, M. C. deAndrade, B. A. Higa, and A. M. Leese de Escobar, Characterization of large two-dimensional YBa2Cu3O7–δ SQUID arrays, Supercond. Sci. Technol. 29, 084003 (2016).
  * Mitchell  _et al._ [2019] E. E. Mitchell, K.-H. Müller, W. E. Purches, S. T. Keenan, C. J. Lewis, and C. P. Foley, Quantum interference effects in 1D parallel high-Tc SQUID arrays with finite inductance, Supercond. Sci. Technol. 32, 124002 (2019).
  * Oppenländer  _et al._ [2000] J. Oppenländer, C. Häussler, and N. Schopohl, Non-$\Phi_{0}$-periodic macroscopic quantum interference in one-dimensional parallel Josephson junction arrays with unconventional grating structure, Phys. Rev. B 63, 024511 (2000).
  * Reinel  _et al._ [1994] D. Reinel, W. Dieterich, T. Wolf, and A. Majhofer, Flux-flow phenomena and current-voltage characteristics of Josephson-junction arrays with inductances, Phys. Rev. B 49, 9118 (1994).
  * Gali Labarias  _et al._ [2022a] M. A. Gali Labarias, K.-H. Müller, and E. E. Mitchell, Modeling the transfer function of two-dimensional SQUID and SQIF arrays with thermal noise, Phys. Rev. Applied 17, 064009 (2022a).
  * Gali Labarias  _et al._ [2022b] M. A. Gali Labarias, K.-H. Müller, and E. E. Mitchell, The role of coupling radius in 1-D parallel and 2-D SQUID and SQIF arrays, IEEE Trans. Appl. Supercond. 32, 1600205 (2022b).
  * Pegrum [2023] C. Pegrum, Modelling high-$T_{c}$ electronics, Supercond. Sci. Technol. 36, 05300 (2023).
  * Berggren and de Escobar [2015] S. Berggren and A. L. de Escobar, Effects of spread in critical currents for series- and parallel-coupled arrays of SQUIDs and Bi-SQUIDs, IEEE Trans. Appl. Supercond. 25, 1600304 (2015).
  * Crété  _et al._ [2021] D. Crété, J. Kermorvant, Y. Lemaître, B. Marcilhac, S. Mesoraca, J. Trastoy, and C. Ulysse, Comparison of magnetic field detectors based on SQUID and/or Josephson junction arrays with HTSc, IEEE Trans. Appl. Supercond. 31, 1601006 (2021).
  * Foley  _et al._ [1999] C. P. Foley, E. E. Mitchell, S. K. H. Lam, B. Sankrithyan, Y. M. Wilson, D. L. Tilbrook, and S. J. Morris, Fabrication and characterisation of YBCO single grain boundary step edge junctions, IEEE Trans. Appl. Supercond. 9(2), 4281 (1999).
  * Mitchell and Foley [2010] E. E. Mitchell and C. P. Foley, YBCO step-edge junctions with high $I_{c}$$R_{n}$, Supercond. Sci. Technol. 23, 065007 (2010).
  * Lam  _et al._ [2014] S. K. Lam, J. Lazar, J. Du, and C. P. Foley, Critical current variation in YBa2Cu3O7-x step-edge junction arrays on MgO substrates, Supercond. Sci. Technol. 27, 055011 (2014).
  * Khapaev  _et al._ [2001] M. M. Khapaev, A. Y. Kidiyarova-Shevchenko, P. Magnelind, and M. Y. Kupriyanov, 3D-MLSI: Software package for inductance calculation in multilayer superconducting integrated circuits, IEEE, Trans. Appl. Supercond. 11, 1090 (2001).
  * Khapaev [2005] M. M. Khapaev, Numerical solution of an integro-differential equation for a sheet current, Differential Equations 41, 1019 (2005).
  * Brandt [2005] E. H. Brandt, Thin superconductors and SQUIDs in perpendicular magnetic field, Phys. Rev. B 72, 024529 (2005).
  * Khapaev and Kupriyanov [2010] M. M. Khapaev and M. Y. Kupriyanov, Sheet current model for inductances extraction and Josephson junctions devices simulation, J. Phys.: Conf. Ser. 248, 012041 (2010).
  * Bishop-Van Horn and Moler [2022] L. Bishop-Van Horn and K. A. Moler, SuperScreen: An open-source package for simulating the magnetic response of two-dimensional superconducting devices, Comp. Phys. Comm. 280, 108464 (2022).
  * Clem and Brandt [2005] J. R. Clem and E. H. Brandt, Response of thin-film SQUIDs to applied fields and vortex fields: Linear SQUIDs, Phys. Rev. B 72, 174511 (2005).
  * Gross  _et al._ [1997] R. Gross, L. Alff, A. Beck, O. M. Froehlich, D. Koelle, and A. Marx, Physics and technology of high temperature superconducting Josephson junctions, IEEE Trans. Appl. Superrcond. 7, 2929 (1997).
  * Keenan  _et al._ [2021] S. Keenan, C. Pegrum, M. Gali Labarias, and E. E. Mitchell, Determining the temperature-dependent London penetration depth in HTS thin films and its effect on SQUID performance, Appl. Phys. Lett. 119, 142601 (2021).
  * Tesche and Clarke [1977] C. D. Tesche and J. Clarke, dc SQUID: Noise and optimization, J. Low Temp. Phys. 29, 301 (1977).
  * Pearl [1964] J. Pearl, Current distribution in superconducting films carrying quantized fluxoids, Appl. Phys. Lett. 5, 65 (1964).
  * Müller  _et al._ [2001] J. Müller, S. Weiss, R. Gross, R. Kleiner, and D. Koelle, Voltage-flux-characteristics of asymmetric dc SQUIDs, IEEE Trans. Appl. Supercond. 11, 912 (2001).


# Refined Motion Compensation with Soft Laser Manipulators using Data-Driven
Surrogate Models

Yongjun Yan13, Qingpeng Ding1, Mingwu Li2, Junyan Yan1 and Shing Shin Cheng13
1Department of Mechanical and Automation Engineering
The Chinese University of Hong Kong, Hong Kong SAR, China
Email<EMAIL_ADDRESS>2Department of Mechanics and Aerospace Engineering
Southern University of Science and Technology, 518055 Shenzhen, China  3Multi-
Scale Medical Robotics Center, Hong Kong SAR, China

###### Abstract

Non-contact laser ablation, a precise thermal technique, simultaneously cuts
and coagulates tissue without the insertion errors associated with rigid
needles. Human organ motions, such as those in the liver, exhibit rhythmic
components influenced by respiratory and cardiac cycles, making effective
laser energy delivery to target lesions while compensating for tumor motion
crucial. This research introduces a data-driven method to derive surrogate
models of a soft manipulator. These low-dimensional models offer computational
efficiency when integrated into the Model Predictive Control (MPC) framework,
while still capturing the manipulator’s dynamics with and without control
input. Spectral Submanifolds (SSM) theory models the manipulator’s autonomous
dynamics, acknowledging its tendency to reach equilibrium when external forces
are removed. Preliminary results show that the MPC controller using the
surrogate model outperforms two other models within the same MPC framework.
The data-driven MPC controller also supports a design-agnostic feature,
allowing the interchangeability of different soft manipulators within the
laser ablation surgery robot system.

##  I Introduction

Laser thermal ablation (LTA) and soft laser manipulators offer enhanced
precision and safety during non-contact procedures [11, 8], eliminating
ablation errors resulting from tissue deformation and needle deflection during
insertion in traditional thermal ablation techniques with a rigid needle [4].
However, precisely tracking tumor motion with a soft laser manipulator is
challenging due to the infinite degrees of freedom (DoFs) in the manipulator’s
structure and the consistently deformable nature [11]. Current analytical
models and neural network approaches have limitations [13], requiring models
that balance accuracy and efficiency for effective motion compensation during
laser ablation.

Data-driven modeling with structural priors provides explainable and
parameterizable models for continuum robots, making them suitable for
integration into mature control theories such as Model Predictive Control
(MPC) to ensure robustness, stability, and invariance-based safety properties
of the soft laser manipulator control system [7, 10]. In soft robotics,
Koopman operator theory has been used to develop linear and nonlinear models
for MPC [2]. While linear Koopman models offer an acceptable computational
footprint but lack accuracy, nonlinear Koopman models are more accurate but
computationally intensive for real-time implementation. Contrary to the
Koopman theory, which lifts nonlinear system dynamics into a high-dimensional
latent state space, recent studies have shown that the dynamics of high-
dimensional systems can be learned on low-dimensional, attracting invariant
manifolds called Spectral Submanifolds (SSM) [3, 9]. Preliminary research
indicates that integrating the SSM model reduction method into an optimal
control framework can result in more accurate tracking performance, as
demonstrated with a toy platform [1].

In this work, we introduce an MPC-based, data-driven control framework for
motion compensation using a soft laser manipulator. By applying SSM theory to
a soft robot in a medical context, we develop a low-dimensional dynamic model
from vibration trajectories with and without cable tension to serve as the
surrogate model of the soft laser manipulator (Fig.1 (a)). A tailored solver
utilizing the sequential quadratic programming (SQP) algorithm addresses the
optimal control problem with the learned surrogate dynamic model and
constraints[12]. Suboptimal and warm-start features of the SQP algorithm are
applied to further reduce the computational burden. Experiments demonstrate
that our proposed controller framework with surrogate models achieves superior
tracking accuracy in periodic trajectories mimicking tumor motion induced by
respiration. Additionally, we validate the design-agnostic feature of the
proposed controller in minimally invasive tumor ablation surgeries.

##  II Mechanisms of Soft Laser Manipulator

Figure 1: Schematic diagram of a soft laser manipulator in minimally invasive
liver tumor ablation surgery. (a) Demonstration of the motion compensation
mechanism with the surrogate models of the soft laser manipulator. (b)
Clinical setup of the minimally invasive tumor ablation surgery with soft
laser manipulator. (c) CAD models of the soft laser manipulator robotics
system. (d) Exposed electro-mechanical subunits of one of the cable actuation
systems. The four cable actuation systems are mounted in pairwise opposition
in the cable actuation modules.

The clinical setup of the minimally invasive tumor ablation surgery using the
soft laser manipulator is shown in Fig.1 (b). The surgical robot equipped with
the laser fiber is inserted through a small incision to deliver non-contact
focused laser ablation to a peripheral liver tumor site without penetrating
the liver tissue [6].

The cable-driven soft laser manipulator is designed with a
$D_{s}=14.0~{}\text{mm}$ outer diameter, a $D_{l}=8.4~{}\text{mm}$ lumen
diameter, and a $L_{f}=94.0~{}\text{mm}$ bending length. The end effector,
shown in Fig.1(c), consists of a soft continuum body, a top cap, and an end
cap. The manipulator is made from materials with significant stiffness
differences: the soft body is vacuum-molded from polyurethane elastomer with a
hardness of 70 Shore A, and the top and end caps are 3D-printed from
photosensitive resin. To enable two degrees of freedom, four 1.2 mm diameter
holes are spaced $90^{\circ}$ apart to create channels for the driving cables.
Thirteen diamond-shaped notches (1.2 mm length, 11.0 mm pitch) along each
channel provide anisotropic stiffness and directional deformation. The
actuation setup for one cable is shown in Fig.1(d). A DC motor-driven lead
screw linear stage pulls and releases the cable along the manipulator’s
length, with a load cell measuring tension at the cable’s proximal end. A 1.0
mm diameter hole at the center of the top cap mounts the laser fiber and
Aurora sensor. A 420 µm diameter ball lens focuses the laser light on the
specimen surface. The end cap supports the continuum body under cable tension,
allowing controlled movement.

The diamond-shaped notches along each cable channel are designed for
controllable bending but also introduce anisotropic stiffness and damping
properties, which are challenging to integrate into current dynamic control
paradigms. Additionally, the slender, deformable structure of the soft laser
manipulator results in twisting and shearing under transverse gravity loads
and cable action, leading to kinematic shifts in tip orientation. Precise
motion compensation in laser ablation therapy requires the manipulator’s
controller to manage a high-frequency, small-diameter 3D periodic trajectory,
necessitating specialized control strategies.

##  III Controller Implementation Pipeline

To illustrate the training process of data-driven surrogate models for the
soft manipulator using low-dimensional SSM theory and to explain how SQP
solves the optimal control problem, the implementation process of the MPC
controller with the surrogate model is introduced in a pipeline style, as
shown below.

In Step 1, snapshots are constructed by allowing the soft manipulator to decay
from 18 different initial points to learn its autonomous dynamics. The
Cartesian positions of two feature points, 1 mm apart along the longitudinal
direction at the tip of the soft manipulator, are measured using an Aurora 6
DoF sensor. A time-delay of $d=1$ is applied to the 18 decaying trajectories
from these two feature points to derive the training snapshot.

In Step 2, Singular Value Decomposition (SVD) is performed, and the leading
$n$ modes of the snapshot are selected to capture the majority of the decay
pattern of the soft manipulator dynamics. Each trajectory $\mathbf{T}_{k}$ in
the observable space is then projected onto the leading $n$ modes to obtain
low-dimensional coordinates via the transformation:
$\mathbf{X}_{k}=\tilde{\mathbf{V}}^{T}_{n}\mathbf{T}_{k}$ [3]. The reduced
dynamics on the manifold are learned via ridge regression,

$\mathbf{R}^{*}=\underset{\mathbf{R}}{\arg\min}\sum_{k}\left\|\dot{\mathbf{T}}_{k}^{d}-\mathbf{R}_{1}\mathbf{X}_{k}-\mathbf{R}_{2:n_{r}}\mathbf{X}_{k}^{2:n_{r}}\right\|_{2}^{2},$
(1)

where all the eigenvalues of the linear state matrix $\mathbf{R}_{1}$ have
negative real parts (these eigenvalues are also associated with the singular
values), ensuring that the SSM attracts trajectories to the equilibrium
points.

In Step 3, the control effect is incorporated into the dynamics equation by
regressing the control matrix $\mathbf{C}_{r}$. A new dataset is collected by
generating a smooth sequence of inputs through periodic actuation at various
control amplitudes and recording the corresponding observed state
trajectories. The mapping and updated dynamics equation with the control input
is reorganized as:

$\dot{\mathbf{x}}=\mathbf{r}(\mathbf{x},\mathbf{u})=\mathbf{R}_{1}\mathbf{x}+\mathbf{R}\mathbf{x}^{2:n_{r}}+\mathbf{C}_{\mathbf{r}}\mathbf{u}.$
(2)

where $\mathbf{x}$ and $\mathbf{u}$ represent the state variables in the low-
dimensional coordinates, and control variables, respectively. Considering the
delay-embedding can not be performed with the current state measures, a
mapping from the observable space to the low-dimensional space is fitted to
facilitate the formulation of the optimal control problem:

$\mathbf{x}=\mathbf{v}_{z}(\mathbf{z})=\tilde{\mathbf{V}}_{0}\mathbf{z}+\tilde{\mathbf{V}}\mathbf{z}^{2:n_{v}}$
(3)

where $\mathbf{z}$ represent the state variables in the observable space.

In Step 4, the synthesized data-driven surrogate model is integrated into an
MPC framework to calculate optimal cable tension for motors, steering the
laser dot to track tumor motion. MPC is chosen because it: 1) proactively
adjusts cable tension by considering look-ahead trajectory information,
compensating for lower-level tracking delays, and 2) naturally handles cable
tension constraints. The optimal control problem is formulated as a nonlinear
programming problem with a multiple-shooting method:

$\min_{x(\cdot|j)u(\cdot|j)}\left\|\delta\mathbf{x}\left(N|j\right)\right\|_{\mathbf{Q}_{f}}^{2}+\sum_{i=0}^{N-1}\|\delta\mathbf{x}(i|j)\|_{\mathbf{Q}}^{2}$
(4)

subject to

$\displaystyle\mathbf{x}(0|j)=\mathbf{v}_{z}(\mathbf{z}(0|j)),$ (5)
$\displaystyle{\mathbf{x}}(i+1|j)=\mathbf{x}(i|j)+\mathbf{r}(\mathbf{x}(i|j),\mathbf{u}(i|j))\Delta
t,i=0:N_{h}-1$
$\displaystyle\mathbf{x}_{ref}(i|j)=\mathbf{v}_{z}(\mathbf{z}_{ref}(i|j)),i=0:N_{h}$
$\displaystyle u(i|j)\in\mathcal{U},i=0:N_{h}.$

Here, $\delta\mathbf{x}(i|j)=\mathbf{x}(i|j)-\mathbf{x}_{ref}(i|j)$ represents
the tracking error in low-dimensional coordinates. Unlike the approach in [1],
which formulates the cost function in observable space, this research
formulates the tracking error in low-dimensional space by mapping the
reference trajectory to the low-dimensional space in advance to avoid the
numerical instability of the SQP solver. A tailored solver based on the SQP
algorithm is employed to solve the nonlinear programming problem, leveraging
suboptimal and warm-start features to reduce computational burden.

##  IV Experimental Evaluation

Figure 2: Experimental platform for the soft laser manipulator.

A testbed (Fig.2) was developed to evaluate the motion compensation
performance of the MPC controller with the surrogate model for the soft laser
manipulator. At the far end of the manipulator is a Stewart platform,
replicating tumor motion induced by heartbeat and respiration.

###  IV-A Control Performance Comparison with Other Models

Two other MPC controllers, namely the Constant Curvature Model-based MPC (CC-
MPC) controller [5] and the Linear Koopman Model-based MPC (LK-MPC) controller
[2], are benchmarked to assess the efficacy of integrating the low-dimensional
surrogate model into the MPC scheme (SSM-MPC).

Figure 3: Reference tracking performance comparison of the SSM-MPC (row 1,
column 1, red), LK-MPC (row 1, column 2, blue), and CC-MPC (row 1, column 3,
yellow). Open-loop prediction accuracy comparison using identical initial
states and control sequences (row 2).

As shown in Fig.3, the SSM-MPC controller with the low-dimensional surrogate
model achieved the best tracking performance, with an overall average error of
0.27 mm, followed by the LK-MPC controller with an average error of 0.36 mm.
The 25% tracking error reduction of the SSM-MPC controller is significant, as
it is achieved using a more accurate model prediction with only a three-
dimensional SSM model compared to the 10-dimensional linear Koopman model with
Fourier basis functions (row 2 in Fig.3). This result aligns with the findings
in [2] that the linear Koopman model cannot capture all the nonlinearity of
the soft manipulator. Therefore, the proposed SSM-MPC controller outperforms
the LK-MPC in tracking accuracy while avoiding the computational burden of
nonlinear Koopman models.

The CC-MPC controller exhibited the worst tracking performance with an average
error of 0.58 mm, nearly twice that of the SSM-MPC controller (Fig.3). This
poor performance is due to the inaccuracy of the constant curvature model,
where the constant curvature assumption fails in the experimental setup of the
soft laser manipulator due to factors such as gravity and anisotropic
properties.

###  IV-B Design Agnostic Feature Evaluation

Figure 4: Design-agnostic feature evaluation of the SSM-MPC controller. (a)
Soft manipulator with a rectangular curve and 70 Shore A hardness white
polyurethane elastomer. (b) Learned submanifolds for the manipulator in (a).
(c) Tracking performance of the SSM-MPC controller with the manipulator in
(a). (d) Soft manipulator with a diamond curve and 90 Shore A hardness black
polyurethane elastomer. (e) Learned submanifolds for the manipulator in (d).
(f) Tracking performance of the SSM-MPC controller with the manipulator in
(d).

Compared with the first-principles models used in CC-MPC, the data-driven
surrogate model that encodes the SSM as structural priors not only offers
superior tracking performance for SSM-MPC but also maintains generalization
capability across different soft manipulator designs. Similar to the parameter
identification process for first-principles models, which involves parameters
with physical interpretations, the surrogate model generalizes by learning
different submanifolds with updated snapshots that reflect specific design
variations. This design-agnostic feature allows the SSM-MPC controller to
function effectively for various soft manipulator designs tailored to
different ablation tasks without requiring tedious parameter tuning.

The design agnostic capability of the SSM-MPC controller is assessed and
depicted in Fig.4. With two distinct designs as depicted in Fig.4 (a) and (d),
the SSM-MPC controller adeptly follows the trajectory, maintaining a tracking
error of 0.40 mm. This performance is achieved by employing tailored SSM
models (shown in Fig.4 (b) and (e)) corresponding to each specific design.

##  V Conclusion

The SSM is an ideal mathematical tool for depicting the attenuating
oscillation pattern of the soft manipulator, as it effectively attracts
trajectories to equilibrium points. The data-driven surrogate model leveraging
SSM theory outperforms the linear Koopman and constant curvature models when
integrated into the MPC framework for motion compensation of the soft
manipulator. By learning the submanifolds of different soft manipulator
designs from corresponding datasets, the data-driven MPC controller supports a
design-agnostic feature, enabling the interchangeability of various soft
manipulators within the laser ablation surgery robot system.

## Acknowledgments

Y. Yan acknowledges the support from the Multi-Scale Medical Robotics Center
that made this work possible.

## References

  * Alora et al. [2023] John Irvin Alora, Mattia Cenedese, Edward Schmerling, George Haller, and Marco Pavone.  Practical deployment of spectral submanifold reduction for optimal control of high-dimensional systems.  _IFAC-PapersOnLine_ , 56(2):4074–4081, 2023.
  * Bruder et al. [2020] Daniel Bruder, Xun Fu, R Brent Gillespie, C David Remy, and Ram Vasudevan.  Data-driven control of soft robots using koopman operator theory.  _IEEE Transactions on Robotics_ , 37(3):948–961, 2020.
  * Cenedese et al. [2022] Mattia Cenedese, Joar Axås, Bastian Bäuerlein, Kerstin Avila, and George Haller.  Data-driven modeling and prediction of non-linearizable dynamics via spectral submanifolds.  _Nature communications_ , 13(1):872, 2022.
  * Crocetti et al. [2010] Laura Crocetti, Thierry de Baere, and Riccardo Lencioni.  Quality improvement guidelines for radiofrequency ablation of liver tumours.  _Cardiovascular and interventional radiology_ , 33:11–17, 2010.
  * Della Santina et al. [2023] Cosimo Della Santina, Christian Duriez, and Daniela Rus.  Model-based control of soft robots: A survey of the state of the art and open challenges.  _IEEE Control Systems Magazine_ , 43(3):30–65, 2023.
  * Fusaglia et al. [2016] Matteo Fusaglia, Hanspeter Hess, Marius Schwalbe, Matthias Peterhans, Pascale Tinguely, Stefan Weber, and Huanxiang Lu.  A clinically applicable laser-based image-guided system for laparoscopic liver procedures.  _International journal of computer assisted radiology and surgery_ , 11:1499–1513, 2016.
  * Hewing et al. [2020] Lukas Hewing, Kim P Wabersich, Marcel Menner, and Melanie N Zeilinger.  Learning-based model predictive control: Toward safe learning in control.  _Annual Review of Control, Robotics, and Autonomous Systems_ , 3:269–296, 2020.
  * Kundrat et al. [2019] Dennis Kundrat, Andreas Schoob, Thomas Piskon, René Grässlin, Patrick J Schuler, Thomas K Hoffmann, Lüder A Kahrs, and Tobias Ortmaier.  Toward assistive technologies for focus adjustment in teleoperated robotic non-contact laser surgery.  _IEEE Transactions on medical robotics and bionics_ , 1(3):145–157, 2019.
  * Li et al. [2022] Mingwu Li, Shobhit Jain, and George Haller.  Nonlinear analysis of forced mechanical systemswith internal resonance using spectral submanifolds, part i: Periodic response and forced response curve.  _Nonlinear Dynamics_ , 110(2):1005–1043, 2022.
  * Markovsky et al. [2023] Ivan Markovsky, Linbin Huang, and Florian Dörfler.  Data-driven control based on the behavioral approach: From theory to applications in power systems.  _IEEE Control Systems Magazine_ , 43(5):28–68, 2023.
  * Runciman et al. [2019] Mark Runciman, Ara Darzi, and George P Mylonas.  Soft robotics in minimally invasive surgery.  _Soft robotics_ , 6(4):423–443, 2019.
  * Yan et al. [2023] Yongjun Yan, Nan Li, Jinlong Hong, Bingzhao Gao, Jia Zhang, Hong Chen, Jing Sun, and Ziyou Song.  Eco-coasting controller using road grade preview: Evaluation and online implementation based on mixed integer model predictive control.  _IEEE Transactions on Vehicular Technology_ , 2023.
  * Yasa et al. [2023] Oncay Yasa, Yasunori Toshimitsu, Mike Y Michelis, Lewis S Jones, Miriam Filippi, Thomas Buchner, and Robert K Katzschmann.  An overview of soft robotics.  _Annual Review of Control, Robotics, and Autonomous Systems_ , 6:1–29, 2023.