# Vortex Dynamics and Dissipation Under High-amplitude Microwave Drive

Mattia Checchin<EMAIL_ADDRESS>A. Grassellino  Fermi National Accelerator
Laboratory, Batavia, Illinois 60510, USA

###### Abstract

In this paper, we describe the vortex dynamics under high-amplitude microwave
drive and its effect on the surface resistance of superconductors. The vortex
surface resistance is calculated with a Montecarlo approach, where the vortex
motion equation is solved for a collection of vortex flux lines each
oscillating within a random pinning landscape. This approach is capable of
providing a detailed description of the microscopic vortex dynamics and in
turn important insights into the microwave field amplitude dependence of the
vortex surface resistance. The numerical simulations are compared against
experimental data of vortex surface resistance at high microwave amplitude
measured by means of bulk niobium superconducting-radio frequency cavities
operating at 1.3 GHz. The good qualitative agreement of simulations and
experiments suggests that the non-linear dependence of the trapped flux
surface resistance with the microwave field amplitude is generated by
progressive microwave depinning and vortex jumps.

##  I Introduction

Upon cooldown below critical temperature in presence of an external magnetic
field, magnetic flux quanta—so-called vortexes—can exist in thermodynamic
equilibrium in the mixed state of type-II superconductors.Abrikosov (1957);
Shubnikov  _et al._ (1937)

Below the lower critical field vortexes are not stable in the superconductor,
however, because of the occurrence of defects in real materials, vortexes get
pinned and survive even in the Meissner state. Incomplete Meissner effect
where 100% of the applied field is trapped by the superconductor can be
achieved with superconducting radio-frequency (SRF) cavities by field-cooling
slowly across transition, as experimentally observed by many studies.Romanenko
_et al._ (2014a, b); Posen  _et al._ (2016); Martinello  _et al._ (2016)

Driven by radio-frequency (RF) and microwave fields, vortexes oscillate
dissipating power and contributing to the total surface resistance. Other than
vortex oscillation, others contributions to the surface resistance are: the
temperature dependent surface resistance due to thermally-excited quasi-
particles above the gap,Mattis and Bardeen (1958) surface resistance due to
sub-gap states,Gurevich (2017) proximity-coupled normal-conducting inclusions
(or oxides),Kubo and Gurevich (2019) nano-hydrides,Romanenko  _et al._ (2013a)
and dielectric dissipation due to spurious two-level systems at low microwave
amplitudes and low temperatures.Romanenko and Schuster (2017); Romanenko  _et
al._ (2018)

Figure 1: In graph (a), the experimental data of trapped flux surface
resistance up to high values of peak magnetic field is reported, while the
normalized $R_{fl}$ over the trapped magnetic field is reported in graph (b).

Vortex dissipation under microwave drive is a topic that has been under
investigation for many decades. Initial studies on type-II superconducting
filmsCardona  _et al._ (1970) brought about the development of models to
describe the vortex surface resistance dependence on frequency.Gittleman and
Rosenblum (1966); Rabinowitz (1971) Since then, several other studies
investigated the microwave behavior of conventional superconductors,Chin  _et
al._ (1992); Janjušević  _et al._ (2006); Alimenti  _et al._ (2019)
cuprates,Revenaz  _et al._ (1994); Matsuda  _et al._ (1995); Golosovsky  _et
al._ (1996) and iron-based superconductorsHashimoto  _et al._ (2009a, b);
Okada  _et al._ (2012) in presence of vortexes. In parallel, many theoretical
works allowed to describe the expected penetration depth change due to vortex
oscillation,Coffey and Clem (1991) the vortex surface impedance temperature-
dependence in high-$\text{T}_{\text{c}}$ superconductors,Marcon  _et al._
(1991) and the precise electrodynamics of vortexes lattices in the mixed
state.Sonin  _et al._ (1992)

In the framework of modern particle accelerators employing SRF cavities as
accelerating devices, dissipation due to vortex motion under microwave drive
is a topic of central importance, inasmuch trapped-fields in the order of
$\sim$1% of the earth magnetic field can increase the surface resistance of
about 1 n$\Omega$,Martinello  _et al._ (2016) a large value for state-of-the-
art cavities capable of achieving surface resistance values as low as 4
n$\Omega$ at 2 K.Grassellino  _et al._ (2013); Martinello  _et al._ (2018)

Early studies on S-band SRF bulk Nb cavitiesPiosczyk  _et al._ (1973)
recognized the effect of vortexes on the surface resistance, identifying the
microwave amplitude dependence of the latter to follow a non-linear behavior.
Later studies showed that up to moderate gradients the microwave field
amplitude dependence of the vortex surface resistance in Nb-sputtered on
copper SRF cavities is linear.Benvenuti  _et al._ (1999) Studies of vortex
dynamics in cavity-grade Nb and theoretical description of the phenomenon
observed in the zero microwave field amplitude approximation were also studied
in detail.Gurevich and Ciovati (2013)

More recently, intensive studies of the trapped flux surface resistance at low
and moderate accelerating gradients as a function of the electron mean-free-
path ($l$) discovered a bell-like trend as a function of $l$,Martinello  _et
al._ (2016); Gonnella  _et al._ (2016) which is well described by the
interplay of pinning- and flux-flow-limited dissipation.Checchin  _et al._
(2017, 2018); Calatroni and Vaglio (2017) The studies also reported an almost
linear dependence of the vortex dissipation with the microwave field amplitude
in the resonator up to moderate values.

The microwave field amplitude dependence can be understood as the occurrence
of non-quadratic terms in the pinning potential.Calatroni and Vaglio (2019)
The linear dependence was shown to be well described by adding a cubic term to
a parabolic pinning potential in the vortex motion equation.Calatroni and
Vaglio (2019) However, this approach is limited since built around the single-
vortex dynamics, while in the material there are several vortexes contributing
to the overall dissipation, each one interacting with its surrounding.

A second approach,Liarte  _et al._ (2018) that overcomes the latter issue,
implements a mean-field method to solve the vortex motion equation for
collective weak pinning. In this case the microwave field amplitude dependence
has two regimes: i) hysteretic losses at low amplitudes where the dependence
is linear, and ii) viscous losses at higher amplitudes where the surface
resistance saturates to a constant value.

In the present work, the high-amplitude microwave field dependence of the
vortex surface resistance is explored by comparing the experimental data
acquired for a high-performing single cell Nb accelerating cavity to the
surface resistance calculated by means of numerical simulations. The model
developed is based on a Montecarlo approach that assumes a random distribution
of pinning potentials in the vortex oscillation plane. The simulation scheme
here presented allows to gather insights into the dynamics of vortexes
oscillation in random pinning potentials such us vortex jumps an microwave
depinning, which are not captured by previous models, but key to correctly
interpret the experimental data.

##  II Experimental Data

Experimental data of vortex surface resistance was estimated from microwave
measurements of a single-cell TESLA-typeAune  _et al._ (2000) bulk niobium
cavity operating at 1.3 GHz.

The cavity surface was prepared with a series of treatments designed to
increase the maximum accelerating gradient achievable by the resonator, so
that to capture the power dissipation fingerprint due to the high-amplitude
microwave behavior of vortexes. The cavity was initially electropolished to
remove 200 $\mu$m and subsequently baked at
$800\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for 3 hours to degas hydrogen and
release internal stress. After this initial surface preparation the cavity was
high-pressure water rinsed in a class 10 clean-room to remove field-emitting
particles form the inner surface and assembled with antennae needed for
microwave characterization. The cavity was then evacuated and baked at
$120\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for 48 hours in situ.

The resonator was equipped with three single-axis flux gates magnetometers
from Bartington Instruments located symmetrically around the cavity equator, a
set of Helmholtz coils aligned coaxially to the cavity beam pipe, and a
temperature mapping system (Tmap)Knobloch  _et al._ (1994) to precisely
monitor the cavity temperature rise during the measurements.

The cavity was initially zero-field-cooled (ZFC) to 1.5 K and data Q-factor
versus accelerating gradient acquired. In such conditions, the resonator
reached maximum peak magnetic field ($B_{p}$) of 204.5 mT (48 MV/m) with
Q-factor of $9.6\cdot 10^{9}$. The same cavity was then field-cooled (FC) very
slowly to 1.5 K four times trapping respectively 7.5 mG, 30 mG, 50 mG, and 100
mG upon the superconducting transition.

As described in detail in Ref. Martinello  _et al._ , 2016, the vortex surface
resistance is calculated by subtracting ZFC data from FC data acquired at the
same temperature. By performing the measurement at 1.5 K, the thermal
contribution to the surface resistance is negligible, and the systematic error
introduced by the accuracy of the temperature control in the dewar minimized.

In Fig. 1(a), the vortex surface resistance ($R_{fl}$) as a function of the
peak magnetic field is reported. The normalized vortex surface resistance over
the magnetic field trapped during the cooldown ($B_{0}$)—so-called trapped-
flux sensitivity—is reported in Fig. 1(b).

Figure 2: A sketch of the cavity section indicating the location of
thermometers in shown in figure (a). Graph (b) shows the ratio between local
magnetic field at the surface over the peak magnetic field at each thermometer
location, while the Tmap showing vortex dissipation collected at $B_{p}=166$
mT is reported in (c).

As shown in Fig. 1(b), once the vortex surface resistance is normalized by
$B_{0}$, the field dependence appears the same within the measurement error
($\sim$10% of $R_{fl}$ value)Melnychuk  _et al._ (2014); Martinello  _et al._
(2016). This means that, within the range of $B_{0}$ and $B_{p}$ investigated,
the peak magnetic field amplitude dependence of $R_{fl}$ is dictated by vortex
dynamics and not by thermal feedback.

In agreement with previous measurement of vortex surface resistance as a
function of accelerating gradient performed in thin film and bulk Nb
cavities,Benvenuti  _et al._ (1999); Martinello  _et al._ (2016) $R_{fl}$
increases almost linearly up to moderate fields ($\sim$80 mT). For higher
field amplitudes the vortex surface resistance deviates from the almost linear
growth and the slope increases with $B_{p}$, leading to higher values of
resistance.

Experimental data was also acquired by means of the Tmap system, composed of
36 boards arranged symmetrically around the cavity, each incorporating 15
thermometers in contact with the outer surface of the cavity as sketched in
Fig. 2(a).

The data was acquired at $B_{p}=166$ mT both in the ZFC experiment and in the
FC experiment with 100 mG applied during the cooldown. The ZFC tmap data was
then subtracted from the FC data in order to obtain the heating pattern due to
trapped-vortexes only. Temperature mapping data is shown in Fig. 2(c). The
board number is reported in the abscissa, while the thermometer number along
the ordinate.

Interestingly, the vortex dissipation shows a specific pattern with heating
concentrated symmetrically with respect to the equator of the cavity (the
equator location corresponds to thermometer number 8), where, instead, the
temperature differential $\Delta T_{out}$ is much lower than neighboring areas
and tends to zero.

In this configuration, the trapped magnetic field $B_{0}$ is parallel to the
cavity beam axis (as shown in Fig. 2(a)) and at the equator the microwave
currents are parallel to the trapped field. Because of the angular dependence
of the Lorentz force that drives vortex oscillation, this geometrical effect
leads to lower dissipation at the equator compared to areas of the cavity
where the flux has an angle almost perpendicular to the microwave currents.

Important to point out that this behavior is intertwined with the modulation
of the microwave currents that varies with the local amplitude of the surface
magnetic field $B$ at the cavity surface. $B$ is almost constant from
thermometer 4 to thermometer 12—surface peak field $B_{p}$ at thermometer 6
and 10, and it decays rapidly to zero from thermometer 3 to 1 and from 13 to
15 (as shown in Fig. 2(b)).

This modulation of the surface field $B$ explains why the vortex dissipation
is almost zero at the cavity iris (thermometers 1 and 15) where flux and
currents are instead close to orthogonality.

##  III Vortex Dynamics

Numerical simulations of vortex dynamics were performed in order to shed light
on the field dependence of the vortex surface resistance.

Let’s assume a superconducting surface at $z=0$ parallel to the plane $xy$,
where the semi-space $z>0$ is filled by superconducting niobium and the semi-
space $z<0$ is vacuum free space. The magnetic field trapped ($B_{0}$) is
parallel to the $z$ axis, and we assume that upon transition the magnetic flux
is subdivided in individual vortexes not interacting with each other
$B_{0}<<B_{c2}$.

The microwave excitation—with frequency $\omega/2\pi$—at the superconducting
surface has magnetic field oscillating along the $x$ axis, inducing current
parallel to the $y$ axis, which decays exponentially in the material along
$z$, with characteristic decay length equal to the penetration depth
$\lambda$.London and London (1935) The amplitude of the microwave field
coincides with the peak magnetic field $B_{p}$.

Neglecting the inertial term, the vortex motion equation has form:

$\eta\dot{u}\left(t,z\right)=\epsilon
u^{\prime\prime}\left(t,z\right)+f_{p}\left(t,z\right)+\gamma
B_{p}\cos\left(\omega t\right)e^{-z/\lambda}\text{,}$ (1)

where $u\left(t,z\right)$ is the vortex displacement as a function of time and
depth $z$, $\dot{u}$ the time partial derivative, $u^{\prime\prime}$ the
second partial derivative along $z$, and $\gamma=\phi_{0}/(\mu_{0}\lambda)$.
$\epsilon$ is the single-vortex line tension,Brandt (1995) defined as:

$\epsilon=\dfrac{\phi_{0}^{2}}{4\pi\mu_{0}\lambda^{2}}\left[\dfrac{1}{2}+\ln\kappa+e^{-0.4-0.8\ln\kappa-0.1\ln^{2}\kappa}\right]\text{,}$
(2)

with $\kappa=\lambda/\xi$. The vortex viscosity $\eta$ is equal to
$\sigma\phi_{0}B_{c2}$,Bardeen and Stephen (1965) and $\sigma$ is the niobium
normal-state conductivity at low temperature. The pinning force per unit of
length is introduced through $f_{p}\left(t,z\right)$.

###  III.1 Pinning Landscape

Vortex pinning in superconductors is explained as the tendency of the system
to minimize the total loss in condensation energy. Being a vortex itself a
singularity in the condensation energy, it interacts with areas in the
material where superconductivity is suppressed or absent in order to minimize
the overall condensation energy loss.

This mechanism intrinsically implies the minimum dimension for a defect to be
an efficient pinning center as the coherence length $\xi$, which represents
the minimum length over which a given change in condensed electron density can
be made.Ginzburg and Landau (1950)

The coherence length in niobium is of the order of $38$ nm or lower when
scattering centers are present. This implies that any defect in the material
with dimension of the order of $\xi$ or larger could pin vortexes. Due to
their small transverse dimension ($\sim 1-5$ nm), grain boundaries and single
dislocations cannot be consider efficient pinning centers, whereas
dislocations tangles, and material non-uniformity can extend for even hundreds
of nanometers, allowing for more efficient pinning.Martinello (2019); Antoine
(2019)

Due to the defective nature of real materials, many defects may pin
simultaneously a single vortex along its length, making virtually impossible
to measure directly the pinning potential, especially for bulk materials.
Theoretical calculationsThuneberg  _et al._ (1984) suggests that isolated
pinning centers can be approximated by Lorentzian functions of the type
$U(r)=-U_{0}/(1+(r/\xi)^{2})$, with $r$ the radial distance form the vortex on
the plane orthogonal to $\phi_{0}\hat{n}$, while experimental measurements on
thin filmsEmbon  _et al._ (2015) showed that real pinning potentials can be
approximated by a summation of several of these Lorentzian functions.

In the case under study, the pinning landscape is defined in the oscillation
plane of vortexes—$xz$, and it is constructed through a random bottom up
approach based on the summation of $N$ Lorentzian pinning potentials described
above:

$\begin{split}&U(x,z)=\\\
&-\sum_{i}^{N}\Bigg{[}\dfrac{U_{i}(b\xi)^{2}}{(b\xi)^{2}+\big{(}x-X_{i}\big{)}^{2}+a_{i}\big{(}z-Z_{i}\big{)}^{2}}\Bigg{]}\text{.}\end{split}$
(3)

Here, $U_{i}$, $X_{i}$, and $Z_{i}$ represent respectively the $i$-th pinning
potential energy and ($x$,$z$) coordinates in the vortex oscillation plane,
the variable $a_{i}$ is a dimensionless parameter that accounts for the
anisotropy of the pinning potential along the $z$ direction, while $b$ is a
multiplicative factor to account for the pinning center dimension. The pinning
force is defined by the summation of the negative derivative along $x$ of the
$i$-th pinning potential, over the total number of pinning centers $N$.

The coordinates $(X_{i},Z_{i})$ and the anisotropy values $a_{i}$ are
generated randomly within a certain domain, usually $-1<X_{i}<1$ $\mu$m,
$0<Z_{i}<3$ $\mu$m, and $0.001<a_{i}<1$. The pinning energy values $U_{i}$ are
generated randomly following a normal distribution centered around an average
value ($U_{ave}$). $U_{ave}$ is defined as the pinning energy that translates
to a maximum pinning force per unit of length ($f_{p,ave}$) comparable to
values found in literature. Usually for niobium the pinning force per unit of
length varies depending on the material impurity content and crystallographic
structure in between $10^{-6}-10^{-4}$ N/m. Allen and Claassen (1989); Park
_et al._ (1992)

The input value $f_{p,ave}$, together with the pinning potential depth
distribution normal deviation ($\sigma_{U}=5\%U_{ave}$), the coordinates
$(X_{i},Z_{i})$, and parameters $a_{i}$ and $b$, describe completely the
random pinning landscape used in the calculation.

###  III.2 Numerical Simulations

Numerical simulations of the vortex motion are defined within the $xz$ plane
domain where the pinning force landscape is defined by means of Eq. 3,
assuming distributions of pinning energy, location, and anisotropy as
described above.

The initial condition and the boundary conditions used in the calculation are
the following:

$\begin{cases}&u(0,z)=0\\\ &u^{\prime}(t,0)=0\\\
&u^{\prime}(t,z_{max})=0\text{,}\end{cases}$ (4)

where $z_{max}$ represent the maximum depth of the solution domain defined
from 0—the microwave surface, to $z_{max}$.

In several works,Gurevich and Ciovati (2013); Liarte  _et al._ (2018) the
boundary condition at the far edge of the solution domain is set to be
$u(t,z_{max})=0$, implying a vortex strongly pinned in the bulk. On the
opposite, in this work vortexes are allowed to find their own equilibrium
position for any randomly-generated pinning landscape and we set
$u^{\prime}(t,z_{max})=0$. As we will see later, this boundary conditions
choice is mandatory to allow the vortex line displacement in the bulk to
evolve with time until an equilibrium position is found, condition reached
after several oscillations when the simulation converges to a steady-state
solution.

Figure 3: Example of vortex dynamics simulation. In (a), the phase portrait of
the vortex tip is reported. The arrows indicate the direction of the vortex
motion. In light blue the transition from initial condition to steady-state is
shown. In (b), the time evolution of the vortex displacement with respect its
initial position—dashed line—is reported. The solid lines show the vortex
motion in the steady-state, while the dotted lines the time evolution from
initial condition to steady-state.

The vortex motion equation is rewritten to be solved numerically, by defining
lengths in units of $\lambda$ and times in unit of $1/f$:
$u=\tilde{u}\lambda$, $z=\tilde{z}\lambda$, $X_{i}=\tilde{X_{i}}\lambda$,
$Z_{i}=\tilde{X_{i}}\lambda$, and $t=\tilde{t}/f$. The motion equation then
becomes:

$\begin{split}\eta\lambda
f\dot{\tilde{u}}\left(\tilde{t},\tilde{z}\right)&=\dfrac{\epsilon}{\lambda}\tilde{u}^{\prime\prime}\left(\tilde{t},\tilde{z}\right)+f_{p}\left(\tilde{t},\tilde{z}\right)+\\\
&+\gamma B_{p}\text{cos}(2\pi\tilde{t})e^{-\tilde{z}}\text{.}\end{split}$ (5)

The motion equation in Eq. 5 was solved numerically by means of the method of
lines. The spatial derivatives along $z$ are discretized into a grid of $M$
points with a finite elements method, while the time variable is kept
continuous. Equation 5 was thus approximated by a system of $M$ differential
equations each describing the time evolution of $\tilde{u}$ evaluated on one
of the $M$ grid points. To solve Eq. 5, $M$ was chosen in between 1000 and
10000 depending on the complexity of the pinning landscape.

The calculation is let evolve for several periods of oscillation, until the
solution converges to a steady-state condition, identified by a close
trajectory in the vortex motion phase portrait, $\dot{\tilde{u}}$ versus
$\tilde{u}$.

In Fig. 3(a), an example of the vortex motion phase portrait at the surface
($z=0$) calculated for surface peak field $B_{p}=50$ mT, $f_{p,ave}=50$
$\mu$N/m, and $l=300$ nm is shown. The motion starts with the vortex in the
initial condition ($\tilde{u}(0,\tilde{z})=0$) and evolves with time until the
steady-state condition is reached.

In Fig. 3(b), the time evolution of the vortex displacement within the pinning
landscape is shown. The vortex initiates its motion from the initial condition
represented by the straight dashed line at $x=0$. Then, it starts oscillating
at the surface while moving towards negative $x$ values attracted by a nearby
pinning center (darker region). Because of this, the initial motion of the
vortex is chaotic—as shown by the phase space in Fig. 3(a) and by the dotted
lines in Fig. 3(b)—and eventually converges to steady-state oscillations when
the vortex find its equilibrium position attracted by a pinning center, as
shown by the solid lines in Fig. 3(b) on the left side of the solution domain.

As briefly described above, the self-stabilization of the vortex to
equilibrium according to the pinning landscape is allowed by the boundary
condition chosen at the far edge of the solution domain. The condition
$\tilde{u}^{\prime}(\tilde{t},\tilde{z}_{max})=0$ allows the vortex to move
and reach an equilibrium position without being constrained to a specific
location.

Simulations were also performed preparing the initial condition by letting the
vortex relaxing to equilibrium before applying the driving force. The
employment of such simulation strategy resulted in an equivalent outcome as if
$\tilde{u}(0,\tilde{z})=0$ was used, but required longer time to converge to a
steady-state solution. Therefore, all the simulations here presented were
performed using the initial condition reported in Eq. 4.

##  IV Vortex Surface Resistance

Figure 4: Example of vortex surface resistance simulation calculated for
$l=300$ nm and $f_{p,ave}=200$ $\mu$N/m as a function of peak field. In graph
(a), the average single-vortex power dissipation ($\langle P^{(1)}\rangle$)
calculated as a function of the peak magnetic field is reported. In (b) and
(c), the average phase portraits calculated over 100 vortex simulations at
0.001 mT and 300 mT respectively is reported. The parameter $u_{0}$ represent
the $x$ offset of the oscillation center. Thin colored solid line represents
the individual vortex phase portrait, while the thick solid black line their
average. In (d), the simulated $R_{fl}/B_{0}$ curve is shown with abscissa in
semi-logarithmic scale. In the inset, the same curve in the range
$1\text{~{}n}\Omega\text{/mG}\leq R_{fl}/B_{0}\leq
1.8\text{~{}n}\Omega\text{/mG}$ is shown with abscissa in linear scale.

During a FC experiment, we assume that all the magnetic flux applied during
cooldown gets trapped in the material and divided into $N_{v}$ vortexes. Once
an electromagnetic field is applied to the material, inducing currents
parallel to the surface, vortexes oscillates absorbing part of the
electromagnetic energy dissipating power.

The power per unit of surface dissipated by vortex motion is defined as:

$\dfrac{dP}{d\Sigma}=\dfrac{1}{2\mu_{0}^{2}}R_{fl}B_{p}^{2}\text{.}$ (6)

Considering a situation in which the magnetic field $B_{p}$ and $R_{fl}$ are
constant over the surface, $R_{fl}$ can be defined as:

$R_{fl}=\dfrac{2\mu_{0}^{2}}{B_{p}^{2}\Sigma}P\text{,}$ (7)

where $P$ represents the total power dissipated by the collection of the
$N_{v}$ oscillating vortexes.

Figure 5: Example of single-vortex dynamics as a function of the microwave
field. In (a), the single-vortex power dissipation as a function of the peak
magnetic field is shown. In (b), (c), and (d), phase portraits and vortex
displacements as a function of time calculated at the field level identified
by red points and labels 1, 2, and 3 respectively in graph (a) are shown.

Because of the random nature of the pinning landscape in real materials, we
assume that each vortex oscillates within its local pinning landscape, which
differs from vortex to vortex except for its average pinning force. Under this
assumption, the power $P$ represents the sum of the single-vortex dissipated
powers $P_{i}^{(1)}$ calculated over the total number of $N_{v}$ vortexes:

$\begin{split}P&=\sum_{i}^{N_{v}}P_{i}^{(1)}\\\ P_{i}^{(1)}&=\dfrac{\gamma
B_{p}}{T}\int_{0}^{T}\int_{0}^{z_{max}}\dot{u}_{i}(t,z)\text{cos}(\omega
t)e^{-z/\lambda}dzdt\text{,}\end{split}$ (8)

with $\dot{u}_{i}(t,z)$ the $i$-th vortex velocity and $T$ the microwave field
oscillation period.

The number of vortexes trapped in the material is estimated by assuming that
each vortex carries a single flux quantum $\phi_{0}$ and that vortexes are
equally distributed over the area $\Sigma$. In turn, the applied magnetic flux
during the cooldown $B_{0}\Sigma$ is divided into $N_{v}=B_{0}\Sigma/\phi_{0}$
vortexes. Thus Eq. 7 can be rewritten as:

$\begin{split}R_{fl}&=\dfrac{2\mu_{0}^{2}B_{0}}{\phi_{0}B_{p}^{2}}\dfrac{\sum_{i}^{N_{v}}P_{i}^{(1)}}{N_{v}}\\\
&=\dfrac{2\mu_{0}^{2}B_{0}}{\phi_{0}B_{p}^{2}}\langle
P^{(1)}\rangle\text{,}\end{split}$ (9)

where $\langle P^{(1)}\rangle$ represents the average of the $N_{v}$ vortexes
dissipated power.

By applying the same definitions used in Eq. 5, we can rewrite the single-
vortex average dissipated power as:

$\begin{split}\langle P^{(1)}\rangle&=\dfrac{\phi_{0}B_{p}\lambda
f}{\mu_{0}N_{v}}\times\\\
&\times\sum_{i}^{N_{v}}\int_{0}^{1}\int_{0}^{\tilde{z}_{max}}\dot{\tilde{u}}_{i}\left(\tilde{t},\tilde{z}\right)\text{cos}\left(2\pi\tilde{t}\right)e^{-\tilde{z}}d\tilde{z}d\tilde{t}\text{.}\end{split}$
(10)

A Montecarlo approach was implemented to estimate $\langle P^{(1)}\rangle$ by
considering a statistical sample of 100 vortexes. In the range of field
probed—tens to hundreds of mG—and areas considered—order of the square meter,
vortexes exceeds the hundred million unities, making Eq. 10 impossible to
compute by means of a desktop computer if the total number of vortexes
involved was considered.

The vortex motion equation in Eq. 5 was computed for each of the 100 vortexes
generating each time an unique random pinning landscape. $P^{(1)}_{i}$ was
calculated once the solution to Eq. 5 approached the steady-state. The
estimation of $\langle P^{(1)}\rangle$ was carried out by averaging the
outcome of the single-vortex simulations. An example of calculation for
$l=300$ nm and $f_{p,ave}=200$ $\mu$N/m is reported in Fig. 4(a).

The colored solid lines represents the single-vortex power $P^{(1)}_{i}$ as a
function of $B_{p}$ calculated each time for a random pinning landscape, for a
total of 100 vortexes. The solid line with black dots represents instead the
average single-vortex dissipated power $\langle P^{(1)}\rangle$.

In Fig. 4(b) and (c) the average phase portrait of vortex oscillation
calculated over 100 vortexes, respectively at $B_{p}=0.001$ mT and 300 mT, is
shown by the thick black solid line. The thin coloured lines represent the
individual single-vortex phase portraits. Fig. 4(d) and its inset show instead
the calculated normalized vortex surface resistance with abscissa in semi-
logarithmic and linear scales, respectively.

At very low microwave fields, the vortex performs oscillations with small
amplitude, since well constrained inside a local pinning potential. In this
regime, the pinning potential can be approximated by a parabolic function, and
the vortex motion equation (Eq. 5) becomes linear. Thus, the calculated vortex
response is linear, implying an elliptical phase portrait and constant surface
resistance (Fig. 4(b) and (d)).

At high microwave fields, the vortex surface resistance tends to saturate to a
constant value equal to the surface resistance in absence of pinning. In this
case, even if each vortex shows a non-linear response, the average phase
portrait tends to an ellipse (Fig. 4(c)), meaning that the global response
becomes linear for higher microwave fields. Pinning is overcame when vortexes
are driven at higher microwave fields and the surface resistance saturates to
a constant value equal to the vortex surface resistance in absence of pinning,
as shown in the inset of Fig. 4(d).

Such saturation behavior of the vortex surface resistance above a certain
microwave field level was also predicted as generated by viscous losses in
Ref. Liarte  _et al._ , 2018.

Figure 6: Comparison of simulated normalized vortex surface resistance with
experimental data. In (a), simulations for different average pinning force at
fixed mean-free-path and fixed factor $b$. In (b), simulations as a function
of the mean-free-path value for fixed average pinning force and factor $b$. In
(c), simulations as a function of the factor $b$ for fixed average pinning
force and mean-free-path.

In summary, for low microwave fields pinning can be treated as a parabolic
potential—the pinning force is linear and proportional to the vortex
displacement, while for high microwave fields pinning becomes negligible. In
both cases, the vortex dynamics is approximated to be linear leading to a
constant surface resistance.

##  V Microwave Depinning

Simulations of single-vortex dynamics allow to shed light into the mechanism
behind the dependence of $R_{fl}$ against the amplitude of the microwave
drive.

In Fig. 5, an example of single-vortex dynamics is reported. The single-vortex
power dissipation ($P^{(1)}_{i}$) curve shown in Fig. 5(a) was selected from
Fig. 4(a). The curve presents a step like behavior (see highlighted points 1
and 2) as a function of $B_{p}$. Such effect is consequence of a “vortex
jump”, where the vortex gets depinned from a location and captured by a
stronger pinning center nearby, as a consequence of microwave depinning.

Simulations in Fig. 5(b), (c), and (d) show the vortex jump behavior just
discussed and the microwave depinning phenomenon. The darker blue line in the
phase portraits represents steady-state oscillations, while the simulations
underneath show the time evolution of the vortex oscillation in the steady-
state only (for clarity). The numbers 1, 2, and 3 indicate the $B_{p}$ value
at which the simulations were done (63 mT, 75 mT, and 170 mT respectively) and
highlight the respective $P_{i}^{(1)}$ values in Fig. 5(a).

At $B_{p}=63$ mT (point 1), the vortex tip oscillates almost depinned from the
local pinning center. As soon as the $B_{p}$ is increased to $75$ mT (point
2), the oscillation widens and the vortex starts to interact with a nearby
pinning center, which eventually captures it.

Comparing the two simulations (point 1 and 2), appears clear that because of
the jump from one pinning center to a stronger one, the amplitude of the
steady-state oscillations are reduced and so it is the oscillation velocity.
As a result, the power dissipation is reduced accordingly, generating this
abrupt step characteristic in the $B_{p}$ dependence of $P^{(1)}_{i}$.

The progressive microwave depinning phenomenon is appreciable by comparing the
simulations of point 2 and 3. At $B_{p}=75$ mT (point 2) the vortex
oscillation is well constrained by the pinning center, while for higher peak
magnetic field values (e.g 170 mT, point 3) the oscillation widens
substantially reaching areas outside the local pinning potential. If the
conditions are right, vortex jump induced by microwave depinning can repeat at
higher fields, as shown in Fig. 5(a) around 180 mT.

Each vortex might experience jumps and depinning at different $B_{p}$ values
depending on the local pinning landscape, producing an unique $P_{i}^{(1)}$
versus $B_{p}$ characteristic. In the global picture, the microwave amplitude
dependence of $R_{fl}$ is defined by the collective vortex motion and the
abrupt steps expected in the $P_{i}^{(1)}$ due to individual vortex jumps are
averaged out.

##  VI Comparison with Experimental Data

###  VI.1 Microwave Data

In Fig. 6(a), (b), and (c), the trapped flux surface resistance measured
experimentally is compared to several simulations of trapped flux surface
resistance calculated for increasing values of average pinning force
$f_{p,ave}$, mean-free-path $l$, and pinning center size $b\xi$, respectively.

As shown in Fig. 6(a), for larger average pinning force, the $B_{p}$ value at
which saturation is reached increases and the overall $R_{fl}$ curve is
shifted to lower values. The effect of different mean-free-path is instead
shown in Fig. 6(b). The flux-flow viscosity decreases for low $l$ values, the
vortex moves more freely and the saturation value increases. The slope of the
$R_{fl}$ microwave field amplitude dependence is instead dictated by the
dimension of the pinning centers, as shown in Fig. 6(c). The smaller the
pinning center, the steeper the trend and the earlier the saturation is
reached.

Since the pinning landscape is unknown a-priori, it is virtually impossible to
generate quantitative simulations of vortex surface resistance as a function
of the microwave field amplitude. Nevertheless, the simulations can
qualitatively describe the experimental data. Based on the simulations
results, progressive microwave depinning is the most likely culprit for the
microwave field dependence of the vortex surface resistance, the higher the
driving force, the less constrained the vortexes and wider the oscillations,
leading to higher resistance.

orHowever, the experimental data does not present saturation at high $B_{p}$
values, contradicting the simulations. This might occur because: i) saturation
is reached at higher fields, and ii) others non-linearity play a role in the
vortex dynamics. The first possibility is related to the fact that pinning
landscape and mean-free-path used in the simulations may not be accurate
enough—pinning landscape is unknown a-priori and cavities treated as the one
under study are characterized by $l$ varying with depth, with much shorter
values at the surface.Romanenko  _et al._ (2014c) The second possibility is
instead substantiated by the fact that in specific conditions the vortex flux-
flow viscosity decreases with the vortex velocity, introducing another non-
linearity in the motion equation. Such effect, first described by Larkin and
Ovchinnikov,Larkin and Ovchinnikov (1975) will be discussed in Section VII.

###  VI.2 Geometrical Considerations

Due to the not planar geometry of the SRF cavity, the vortex surface
resistance is not constant across the surface, it instead varies as described
previously in Section II.

The temperature variation $\Delta T_{out}$—calculated as the FC and ZFC Tmap
data subtraction—is directly related to the surface resistance due to
vortexes. The surface resistance is directly proportional to the power
dissipated, as shown in Eq. 7, which, in turn, is directly proportional to the
temperature variation on the outside surface of the cavity through the thermal
transport law, i.e. $R_{fl}\propto\Delta T_{out}/B_{p}^{2}$.

Figure 7 reports the normalized data shown in Fig. 2(c) averaged over the
board number as a function of the angle $\theta$ between $B_{0}\hat{z}$ and
the microwave currents flowing parallel to the cavity surface. $\theta$ was
estimated as the angle between $B_{0}\hat{z}$ and the cavity surface tangent
at each thermometer location.

The error bars are defined to be the standard deviation of the $\Delta
T_{out}/B_{p}^{2}$ average value in ordinate, while for the values in abscissa
the error was estimated a-posteriori to be $10\degree$.

The experimental data is compared against the model here presented, where the
Lorentz force $\text{sin}(\theta)$ dependence is taken into consideration, as
well as the local surface current distribution at the cavity surface reported
in Fig. 2(b).

Figure 7: Normalized experimental data $\Delta T/B_{p}^{2}$ and normalized
simulated $R_{fl}$ plotted against the angle between microwave currents and
magnetic flux directions.

The calculation was performed with the same Montecarlo approach described
above, where every point is calculated averaging the dissipated power over 100
vortexes. The simulation was carried out assuming $l=100$ nm, $f_{p,ave}=200$
$\mu$N/m, and $b=4$, but any other set of parameters returned comparable
results. In order to compare experimental data and simulation, all the
quantities in ordinate are normalized between 0 and 1.

The angle dependence calculated is in good agreement to the experimental data.
At $\theta=0$°, the surface resistance is zero since the trapped field is
parallel to the microwave currents, while at $\theta\simeq 80$°, the surface
resistance decreases because the microwave current drops rapidly to zero
approaching the cavity irises.

##  VII Larkin-Ovchinnikov Instability

Larkin-Ovchinnikov (LO) instabilityLarkin and Ovchinnikov (1975) is a
phenomenon of non-equilibrium related to vortex motion. Inside the vortex core
superconductivity is suppressed and quasi-particles can occupy virtually any
energy level. During vortex motion the energy of quasi-particles increases and
when reaches $\sim\Delta$, quasi-particles diffuse out of the vortex into the
superconductor, the vortex core gets depleted, shrinks in dimension, and the
flux-flow viscosity decreases.

In the LO theory framework, the critical velocity is defined asLarkin and
Ovchinnikov (1975); Grimaldi  _et al._ (2010); Leo  _et al._ (2011)

$v^{*}\propto\tau_{\varepsilon}^{-1/2}\left(1-\dfrac{T}{T_{c}}\right)^{1/4}\text{,}$
(11)

where $\tau_{\varepsilon}^{-1}$ is the rate at which quasi-particles relax
back into the vortex core. The shorter $\tau_{\varepsilon}$, the faster the
vortex must travel to meet the instability onset.

Experimental studies on Nb thin filmsGrimaldi  _et al._ (2010) showed that the
LO instability critical velocity $v^{*}$ is function of the magnetic field
trapped in the sample. The low $B_{0}$ region ($B_{0}\leq 100$ mG) which is of
our interest is still poorly explored. However, it was shown that, in Nb,
$v^{*}$ increases for very low fields,Grimaldi  _et al._ (2010); Dobrovolskiy
_et al._ (2017) reaching values of the order of couple of km/h for $B_{0}\sim
1$ G.

The theory developed by Larkin and Ovchinnikov is strictly applicable to
$T\simeq T_{c}$. However, the underlying decreasing of the flux-flow viscosity
as a function of the vortex velocity is a feature observed also at low
temperatures,Leo  _et al._ (2011) and therefore relevant in Nb SRF cavities
operating at 2 K, especially at high gradients. In order to expand the
applicability of such formulation to low temperatures, the quasi-particle
relaxation rate $\tau_{\varepsilon}^{-1}$ shall be defined as the sum
$\tau_{s}^{-1}+\tau_{r}^{-1}$, respectively scattering and recombination
rates.Leo  _et al._ (2011)

Generally speaking, at low temperatures ($T<T_{c}/2$), $\tau_{s}^{-1}$ and
$\tau_{r}^{-1}$ decrease due to the lower number of phonons and quasi-
particles, respectively.Kaplan  _et al._ (1976) However, if $\tau_{s}^{-1}$ is
governed by scattering with impurities or lattice defects, it becomes
temperature-independent. The occurrence of LO instability should then be more
relevant in clean niobium, where the quasi-particle relaxation time is not
limited by scattering with impurities and results longer.

In the present study, we discuss trapped flux surface resistance in cavities
with low mean-free-path, prepared with a mild baking treatment (Section II),
that was proven to increase the concentration of interstitial oxygen and
vacancies in the material.Ciovati (2004); Romanenko  _et al._ (2013b, 2016) We
can then assume that the relaxation rate is most likely dominated by
scattering with impurities, and therefore that $v^{*}$ increases with $T$ as
expected by Eq. 11.

Hence, if occurring, LO instability would be appreciable at higher microwave
fields, introducing deviations from the expected plateau in the surface
resistance predicted by this work. However, since no saturation in the
experimental data was observed at high microwave field amplitudes, LO
instability cannot be totally ruled out, inasmuch saturation might be masked
by such phenomenon.

##  VIII Conclusions

In this work, the vortex surface resistance as a function of the microwave
field amplitude is simulated by means of a Montecarlo approach aimed to solve
the motion equation for a collection of vortexes each oscillating within a
randomly defined pinning landscape.

Simulations suggest that at low microwave amplitudes vortexes are well
constrained within a local pinning potential and their response is linear,
leading to a constant surface resistance. For increasing field amplitudes
vortexes may oscillate out of their local pinning potential and possibly jump
from their local to a neighbouring pinning center introducing non-linearity.
The effect of such phenomena leads to an increase of the vortex surface
resistance until the average response approaches a linear behavior for higher
field amplitudes, where the average vortex motion becomes insensitive to the
pinning landscape. When this happens the surface resistance saturates to a
constant value.

Slope, saturation value, and saturation onset are closely related to the
pinning landscape—not known a-priori—within which vortexes oscillate.
Quantitative description of $R_{fl}$ as a function of $B_{p}$ is then
virtually impossible. Nevertheless, the experimental data are in qualitative
agreement with the simulations especially the data as a function of the angle
between microwave current and vortex magnetic flux direction.

The experimental data as a function of the microwave field amplitude do not
show saturation as instead expected from simulations. As described above, this
might be due to the not enough accurate choice of pinning parameters and mean-
free-path value in the simulations, or to the onset of LO instability at high
fields.

Concluding, this work allowed to gain insights into vortex dynamics at high
amplitude microwave drive and on the transition from low to high field
amplitudes. The microscopic nature the microwave field dependence of the
vortex surface resistance was presented and attributed to microwave depinning
of vortexes. Appreciable contributions from thermal feedback effects for the
magnetic field values and microwave amplitudes explored
experimentally—$B_{0}<100$ mG and $B_{p}<200$ mT—were instead ruled out.

##  IX Acknowledgements

We thank Alex Romanenko and Martina Martinello for the insightful discussions
that help the development of this work. This work was supported by the United
States Department of Energy, Office of High Energy Physics. Fermilab is
operated by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359
with the United States Department of Energy.

##  Appendix A Penetration Depth and Coherence Length Calculation

Since the numeric calculation is performed for finite vales of electron mean-
free-path ($l$) value, the penetration depth of niobium was calculated by
means of the spatial Fourier transform of the Pippard’s Kernel Pippard (1953)
$K(q,l,T)$, as described by Tinkham, Tinkham (2004) that closely approximates
the BCS Bardeen  _et al._ (1957) results.

For isotropic superconductors, such as niobium, $K(q,l,T)$ is defined as:

$\begin{split}&K(q,l,T)=\dfrac{\xi(l,T)}{\xi_{0}(T)\lambda_{L}(T)^{2}}\bigg{\\{}\dfrac{3}{2(q\xi(l,T))^{3}}+\\\
&+\left[(1+q^{2}\xi(l,T)^{2})\text{tan}^{-1}q\xi(l,T)-q\xi(l,T)\right]\bigg{\\}}\text{,}\end{split}$
(12)

where $\lambda_{L}(T)$ is approximated by the two-fluid model definition
$\lambda_{L}(T)=\lambda_{L}(0)\left[1-(T/T_{c})^{4}\right]^{-1/2}$, Gorter and
Casimir (1934) and the $T$ dependence of $\xi$ and $\xi_{0}$ is derived from
the BCS theory: Bardeen  _et al._ (1957)

$\begin{split}\xi(l,T)&=\left[\dfrac{J(0,T)}{\xi_{0}(0)}+\dfrac{1}{l}\right]^{-1}\\\
\xi_{0}(T)&=\dfrac{\xi_{0}(0)}{J(0,T)}\text{.}\end{split}$ (13)

Here, $J(0,T)$ is the real-space Kernel as defined as: Bardeen  _et al._
(1957); Suter  _et al._ (2005)

$J(0,T)=\left[\dfrac{\lambda_{L}(T)}{\lambda_{L}(0)}\right]^{2}\dfrac{\Delta(T)}{\Delta(0)}\text{tanh}\dfrac{\Delta(T)}{2\kappa_{B}T}\text{,}$
(14)

where
$\Delta(T)/\Delta(0)=\left[\text{cos}\left(\pi/2\left(T/T_{c}\right)^{2}\right)\right]^{1/2}$,Sheehan
(1966) and $\Delta(0)$ the energy gap at $T=0$ K.

The penetration depth $\lambda(T)$ is then calculated numerically for specular
reflection at the boundary as:

$\lambda(l,T)=\dfrac{2}{\pi}\int_{0}^{\infty}\dfrac{dq}{K(q,l,T)+q^{2}}\text{.}$
(15)

The parameters $\lambda_{L}(0)$ and $\xi_{0}(0)$ represent the London
penetration depth and the coherence length in the clean limit for $T=0$,
respectively 39 nm and 38 nm.Maxfield and McLean (1965)

## References

  * Abrikosov (1957) A. A. Abrikosov, Zh. Eksp. Teor. Fiz. 32, 1442 (1957), [A. A. Abrikosov, Soviet Phys.—JETP 5, 1174 (1957)].
  * Shubnikov  _et al._ (1937) L. V. Shubnikov, V. I. Khotkevich, Y. D. Shepelev, and Y. N. Riabinin, Zh. Eksp. Teor. Fiz. 7, 221 (1937).
  * Romanenko  _et al._ (2014a) A. Romanenko, A. Grassellino, O. Melnychuk, and D. A. Sergatskov, J. Appl. Phys. 115, 184903 (2014a).
  * Romanenko  _et al._ (2014b) A. Romanenko, A. Grassellino, A. C. Crawford, D. A. Sergatskov, and O. Melnychuk, Appl. Phys. Lett. 105, 234103 (2014b).
  * Posen  _et al._ (2016) S. Posen, M. Checchin, A. C. Crawford, A. Grassellino, M. Martinello, O. S. Melnychuk, A. Romanenko, D. A. Sergatskov, and Y. Trenikhina, J. Appl. Phys 119, 213903 (2016).
  * Martinello  _et al._ (2016) M. Martinello, A. Grassellino, M. Checchin, A. Romanenko, O. Melnychuk, D. A. Sergatskov, S. Posen, and J. F. Zasadzinski, Appl. Phys. Lett. 109, 062601 (2016).
  * Mattis and Bardeen (1958) D. C. Mattis and J. Bardeen, Phys. Rev. 111, 412 (1958).
  * Gurevich (2017) A. Gurevich, Supercond. Sci. Technol. 30, 034004 (2017).
  * Kubo and Gurevich (2019) T. Kubo and A. Gurevich, Phys. Rev. B 100, 064522 (2019).
  * Romanenko  _et al._ (2013a) A. Romanenko, F. Barkov, L. D. Cooley, and A. Grassellino, Supercond. Sci. Technol. 26, 035003 (2013a).
  * Romanenko and Schuster (2017) A. Romanenko and D. Schuster, Phys. Rev. Lett. 119, 264801 (2017).
  * Romanenko  _et al._ (2018) A. Romanenko, R. Pilipenko, S. Zorzetti, D. Frolov, M. Awida, S. Belomestnykh, S. Posen, and A. Grassellino, (2018), arXiv:1810.03703 .
  * Cardona  _et al._ (1970) M. Cardona, G. Fischer, and B. Rosenblum, Phys. Rev. Lett. 12, 101 (1970).
  * Gittleman and Rosenblum (1966) J. I. Gittleman and B. Rosenblum, Phys. Rev. Lett. 16, 734 (1966).
  * Rabinowitz (1971) M. Rabinowitz, J. Appl. Phys. 42, 88 (1971).
  * Chin  _et al._ (1992) C. C. Chin, D. E. Oates, G. Dresselhaus, and M. S. Dresselhaus, Phys. Rev. B 45, 4788 (1992).
  * Janjušević  _et al._ (2006) D. Janjušević, M. S. Grbić, M. Požek, A. Dulčić, D. Paar, B. Nebendahl, and T. Wagner, Phys. Rev. B 74, 104501 (2006).
  * Alimenti  _et al._ (2019) A. Alimenti, N. Pompeo, K. Torokhtii, T. Spina, R. Flükiger, L. Muzzi, and E. Silva, IEEE Trans. Appl. Supercond. 29, 3500104 (2019).
  * Revenaz  _et al._ (1994) S. Revenaz, D. L.-L. D. E. Oates, G. Dresselhaus, and M. S. Dresselhaus, Phys. Rev. B 50, 1178 (1994).
  * Matsuda  _et al._ (1995) Y. Matsuda, M. B. Gaifullin, K. Kumagai, K. Kadowaki, and T. Mochiku, Phys. Rev. Lett. 75, 4512 (1995).
  * Golosovsky  _et al._ (1996) M. Golosovsky, M. Tsindlekht, and D. Davidov, Supercond. Sci. Technol. 9, 1 (1996).
  * Hashimoto  _et al._ (2009a) K. Hashimoto, T. Shibauchi, T. Kato, K. Ikada, R. Okazaki, H. Shishido, M. Ishikado, H. Kito, A. Iyo, H. Eisaki, S. Shamoto, and Y. Matsuda, Phys. Rev. Lett. 102, 017002 (2009a).
  * Hashimoto  _et al._ (2009b) K. Hashimoto, T. Shibauchi, S. Kasahara, K. Ikada, S. Tonegawa, T. Kato, R. Okazaki, C. J. van der Beek, M. Konczykowski, H. Takeya, K. Hirata, T. Terashima, and Y. Matsuda, Phys. Rev. Lett. 102, 207001 (2009b).
  * Okada  _et al._ (2012) T. Okada, H. Takahashi, Y. Imai, K. Kitagawa, K. Matsubayashi, Y. Uwatoko, and A. Maeda, Phys. Rev. B 86, 064516 (2012).
  * Coffey and Clem (1991) M. W. Coffey and J. R. Clem, Phys. Rev. Lett. 67, 386 (1991).
  * Marcon  _et al._ (1991) R. Marcon, R. Fastampa, M. Giura, and E. Silva, Phys. Rev. B 43, 2940 (1991).
  * Sonin  _et al._ (1992) E. B. Sonin, A. K. Tagantsev, and K. B. Traito, Phys. Rev. B 43, 5830 (1992).
  * Grassellino  _et al._ (2013) A. Grassellino, A. Romanenko, D. A. Sergatskov, O. Melnychuk, Y. Trenikhina, A. C. Crawford, A. Rowe, M. Wong, T. Khabiboulline, and F. Barkov, Supercond. Sci. Tech. 26, 102001 (2013).
  * Martinello  _et al._ (2018) M. Martinello, M. Checchin, A. Romanenko, A. Grassellino, S. Aderhold, S. K. Chandrasekeran, O. Melnychuk, S. Posen, and D. A. Sergatskov, Phys. Rev. Lett. 121, 224801 (2018).
  * Piosczyk  _et al._ (1973) B. Piosczyk, P. Kneisel, O. Stoltz, and J. Halbritter, IEEE Trans. Nucl. Sci. 20, 108 (1973).
  * Benvenuti  _et al._ (1999) C. Benvenuti, S. Calatroni, I. E. Campisi, P. Darriulat, M. A. Peck, R. Russo, and A.-M. Valente, Physica C 316, 153 (1999).
  * Gurevich and Ciovati (2013) A. Gurevich and G. Ciovati, Phys. Rev. B 87, 054502 (2013).
  * Gonnella  _et al._ (2016) D. Gonnella, J. Kaufman, and M. Liepe, J. Appl. Phys 119, 073904 (2016).
  * Checchin  _et al._ (2017) M. Checchin, M. Martinello, A. Grassellino, A. Romanenko, and J. F. Zasadzinski, Supercond. Sci. Technol. 30, 034003 (2017).
  * Checchin  _et al._ (2018) M. Checchin, M. Martinello, A. Grassellino, S. Aderhold, S. K. Chandrasekaran, O. S. Melnychuk, S. Posen, A. Romanenko, and D. A. Sergatskov, Appl. Phys. Lett. 112, 072601 (2018).
  * Calatroni and Vaglio (2017) S. Calatroni and R. Vaglio, IEEE Trans. Appl. Supercond. 27, 3500506 (2017).
  * Calatroni and Vaglio (2019) S. Calatroni and R. Vaglio, Phys. Rev. Accel. Beams 22, 022001 (2019).
  * Liarte  _et al._ (2018) D. B. Liarte, D. Hall, P. N. Koufalis, A. Miyazaki, A. Senanian, M. Liepe, and J. P. Sethna, Phys. Rev. Applied 10, 054057 (2018).
  * Aune  _et al._ (2000) B. Aune, R. Bandelmann, D. Bloess, B. Bonin, A. Bosotti, M. Champion, C. Crawford, G. Deppe, B. Dwersteg, D. Edwards, H. T. Edwards, M. Ferrario, M. Fouaidy, P.-D. Gall, A. Gamp, A. Gössel, J. Graber, D. Hubert, M. Hüning, M. Juillard, T. Junquera, H. Kaiser, G. Kreps, M. Kuchnir, R. Lange, M. Leenen, M. Liepe, L. Lilje, A. Matheisen, W.-D. Möller, A. Mosnier, H. Padamsee, C. Pagani, M. Pekeler, H.-B. Peters, O. Peters, D. Proch, K. Rehlich, D. Reschke, H. Safa, T. Schilcher, P. Schmüser, J. Sekutowicz, S. Simrock, W. Singer, M. Tigner, D. Trines, K. Twarowski, G. Weichert, J. Weisend, J. Wojtkiewicz, S. Wolff, and K. Zapfe, Phys. Rev. ST Accel. Beams 3, 092001 (2000).
  * Knobloch  _et al._ (1994) J. Knobloch, H. Muller, and H. Padamsee, Rev. Sci. Instrum. 65, 3521 (1994).
  * Melnychuk  _et al._ (2014) O. Melnychuk, A. Grassellino, and A. Romanenko, Rev. Sci. Instrum. 85, 124705 (2014).
  * London and London (1935) F. London and H. London, Proc. R. Soc. Lond. A 149, 71 (1935).
  * Brandt (1995) E. H. Brandt, Rep. Prog. Phys. 58, 1465 (1995).
  * Bardeen and Stephen (1965) J. Bardeen and M. J. Stephen, Phys. Rev. 140, A1197 (1965).
  * Ginzburg and Landau (1950) V. L. Ginzburg and L. D. Landau, Zh. Eksp. Teor. Fiz. 20, 1064 (1950), [L. D. Landau, _Collected papers_ (Pergamon Press, New York, 1965)].
  * Martinello (2019) M. Martinello, “Microscopic investigation of flux trapping sites in bulk nb,”  (2019), 19th International Conference on RF Superconductivity.
  * Antoine (2019) C. Z. Antoine, Phys. Rev. Accel. Beams 22, 034801 (2019).
  * Thuneberg  _et al._ (1984) E. V. Thuneberg, J.Kurkijärvi, and D. Rainer, Phys. Rev. B 29, 3913 (1984).
  * Embon  _et al._ (2015) L. Embon, Y. Anahory, A. Suhov, D. Halbertal, J. Cuppens, A. Yakovenko, A. Uri, Y. Myasoedov, M. L. Rappaport, M. E. Huber, A. Gurevich, and E. Zeldov, Sci. Rep. 5, 7598 (2015).
  * Allen and Claassen (1989) L. H. Allen and J. H. Claassen, Phys. Rev. B 39, 2054 (1989).
  * Park  _et al._ (1992) G. S. Park, C. E. Cunningham, B. Cabrera, and M. E. Huber, Phys. Rev. Lett. 68, 1920 (1992).
  * Romanenko  _et al._ (2014c) A. Romanenko, A. Grassellino, F. Barkov, A. Suter, Z. Salman, and T. Prokscha, Appl. Phys. Lett. 104, 072601 (2014c).
  * Larkin and Ovchinnikov (1975) A. I. Larkin and Y. N. Ovchinnikov, Zh. Eksp. Teor. Fiz. 68, 1915 (1975), [A. I. Larkin and Yu. N. Ovchinnikov, _Sov. Phys.-JETP_ , 41, 690 (1976)].
  * Grimaldi  _et al._ (2010) G. Grimaldi, A. Leo, D. Zola, A. Nigro, S. Pace, F. Laviano, and E. Mezzetti, Phys. rev. B 82, 024512 (2010).
  * Leo  _et al._ (2011) A. Leo, G. Grimaldi, R. Citro, A. Nigro, S. Pace, and R. P. Hubener, Phys. rev. B 84, 014536 (2011).
  * Dobrovolskiy  _et al._ (2017) O. V. Dobrovolskiy, V. A. Shklovskij, M. Hanefeld, M. Zörb, L. Köhs, and M. Huth, Supercond. Sci. Technol. 30, 085002 (2017).
  * Kaplan  _et al._ (1976) S. B. Kaplan, C. C. Chi, D. N. Langenberg, J. J. Chang, S. Jafarey, and D. J. Scalapino, Phys. Rev. B 14, 4854 (1976).
  * Ciovati (2004) G. Ciovati, J. Appl. Phys. 96, 1591 (2004).
  * Romanenko  _et al._ (2013b) A. Romanenko, C. J. Edwardson, P. G. Coleman, and P. J. Simpson, App. Phys. Lett. 102, 232601 (2013b).
  * Romanenko  _et al._ (2016) A. Romanenko, Y. Trenikhina, M. Martinello, D. Bafia, and A. Grassellino, in  _Proc. of the 19th International Conference on RF Superconductivity_ , THP014 (2016) p. 866.
  * Pippard (1953) A. B. Pippard, Proc. R. Soc. Lond. A 216, 547 (1953).
  * Tinkham (2004) M. Tinkham,  _Introduction to Superconductivity_ (Dover Publications, 2004).
  * Bardeen  _et al._ (1957) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957).
  * Gorter and Casimir (1934) C. J. Gorter and H. Casimir, Physica 1, 306 (1934).
  * Suter  _et al._ (2005) A. Suter, E. Morenzoni, N. Garifianov, R. Khasanov, E. Kirk, H. Luetkens, T. Prokscha, and M. Horisberger, Phys. Rev. B 72, 024506 (2005).
  * Sheehan (1966) T. P. Sheehan, Phys. Rev. 149, 368 (1966).
  * Maxfield and McLean (1965) B. W. Maxfield and W. L. McLean, Phys. Rev. 139, A1515 (1965).


11institutetext: Chair of Process and Data Science, RWTH Aachen University,
Aachen, Germany

# Towards Quantifying Privacy in Process Mining

Majid Rafiei 0000-0001-7161-6927 Wil M.P. van der Aalst  0000-0002-0955-6940

###### Abstract

Process mining employs event logs to provide insights into the actual
processes. Event logs are recorded by information systems and contain valuable
information helping organizations to improve their processes. However, these
data also include highly sensitive private information which is a major
concern when applying process mining. Therefore, privacy preservation in
process mining is growing in importance, and new techniques are being
introduced. The effectiveness of the proposed privacy preservation techniques
needs to be evaluated. It is important to measure both sensitive data
protection and data utility preservation. In this paper, we propose an
approach to quantify the effectiveness of privacy preservation techniques. We
introduce two measures for quantifying disclosure risks to evaluate the
sensitive data protection aspect. Moreover, a measure is proposed to quantify
data utility preservation for the main process mining activities. The proposed
measures have been tested using various real-life event logs.

###### Keywords:

Responsible process mining Privacy preservation Privacy quantification Data
utility Event logs

##  1 Introduction

Process mining bridges the gap between traditional model-based process
analysis (e.g., simulation), and data-centric analysis (e.g., data mining)
[1]. The three basic types of process mining are process discovery, where the
aim is to discover a process model capturing the behavior seen in an event
log, conformance checking, where the aim is to find commonalities and
discrepancies between a process model and an event log, and process re-
engineering (enhancement), where the idea is to extend or improve a process
model using event logs.

An event log is a collection of events. Each event has the following mandatory
attributes: a case identifier, an activity name, a timestamp, and optional
attributes such as resources or costs. In the human-centered processes, case
identifiers refer to individuals. For example, in a patient treatment process,
the case identifiers refer to the patients whose data are recorded. Moreover,
other attributes may also refer to individuals, e.g., resources often refer to
persons performing activities. When event logs explicitly or implicitly
include personal data, privacy concerns arise which should be taken into
account w.r.t. regulations such as the European General Data Protection
Regulation (GDPR).

The privacy and confidentiality issues in process mining are recently
receiving more attention and various techniques have been proposed to protect
sensitive data. Privacy preservation techniques often apply anonymization
operations to modify the data in order to fulfill desired privacy
requirements, yet, at the same time, they are supposed to preserve data
utility. To evaluate the effectiveness of these techniques, their effects on
sensitive data protection and data utility preservation need to be measured.
In principle, privacy preservation techniques always deal with a trade-off
between data utility and data protection, and they are supposed to balance
these aims.

Figure 1: Overview of privacy-related activities in process mining. Privacy
preservation techniques are applied to event logs to provide desired privacy
requirements. The aim is to protect sensitive personal data, yet, at the same
time, preserve data utility, and generate as similar as possible results to
the original ones. The parts indicated by dashed callouts show the analyses
that need to be performed to evaluate the effectiveness of privacy
preservation techniques.

Figure 1 shows the general view of privacy in process mining including two
main activities: Privacy-Preserving Data Publishing (PPDP) and Privacy-
Preserving Process Mining (PPPM). PPDP aims to hide the identity and the
sensitive data of record owners in event logs to protect their privacy. PPPM
aims to extend traditional process mining algorithms to work with the non-
standard event data so-called Event Log Abstraction (ELA) [16] that might
result from PPDP techniques. Abstractions are intermediate results, e.g., a
directly follows graph could be an intermediate result of a process discovery
algorithm. Note that PPPM algorithms are tightly coupled with the
corresponding PPDP techniques.

In this paper, our main focus is on the analyses indicated by the check-boxes
in Fig. 1. Note that disclosure risk analysis is done for a single event log,
while for data/result utility analysis, the original event log/result need to
be compared with the privacy-aware event log/result. We consider simple event
logs containing basic information for performing two main process mining
activities: process discovery and conformance checking. We introduce two
measures for quantifying disclosure risks in a simple event log: identity
(case) disclosure and attribute (trace) disclosure. Using these measures, we
show that even simple event logs could disclose sensitive information. We also
propose a measure for quantifying data utility which is based on the earth
mover’s distance. So far, the proposed privacy preservation techniques in
process mining use the result utility approach to demonstrate the utility
preservation aspect which is not as precise and general as the data utility
approach, since it is highly dependent on the underlying algorithms. We
advocate the proposed measures by assessing their functionality for
quantifying the disclosure risks and data utility on real-life event logs
before and after applying a privacy preservation technique with different
parameters.

The remainder of the paper is organized as follows. Section 2 outlines related
work. In Section 3, formal models for event logs are presented. We explain the
measures in Section 4. The experiments are described in Section 5, and Section
6 concludes the paper.

##  2 Related Work

In process mining, the research field of confidentiality and privacy is
growing in importance. In [2], Responsible Process Mining (RPM) is introduced
as the sub-discipline focusing on possible negative side-effects of applying
process mining. In [12], the authors propose a privacy-preserving system
design for process mining, where a user-centered view is considered to track
personal data. In [18], a framework is introduced providing a generic scheme
for confidentiality in process mining. In [14], the authors introduce a
privacy-preserving method for discovering roles from event data. In [6], the
authors apply $k$-anonymity and $t$-closeness on event data to preserve the
privacy of resources. In [11], the notion of differential privacy is employed
to preserve the privacy of cases. In [17], the $\mathrm{TLKC}$-privacy model
is introduced to deal with high variability issues in event logs for applying
group-based anonymization techniques. In [5], a secure multi-party computation
solution is proposed for preserving privacy in an inter-organizational
setting. In [13], the authors analyze data privacy and utility requirements
for healthcare event data, and the suitability of privacy-preserving
techniques is assessed. In [16], privacy metadata in process mining are
discussed and a privacy extension for the XES standard (https://xes-
standard.org/) is proposed.

Most related to our work is [22], where a uniqueness-based measure is proposed
to evaluate the re-identification risk of event logs. Privacy quantification
in data mining is a well-developed field where the effectiveness of privacy
preservation techniques is evaluated from different aspects such as
dissimilarity [3], information loss [7], discernibility [8], and etc. We
utilize the experiences achieved in this field and propose a trade-off
approach as suggested in [4].

##  3 Preliminaries

In this section, we provide formal definitions for event logs used in the
remainder. An event log is a collection of events, composed of different
attributes, such that they are uniquely identifiable. In this paper, we
consider only the mandatory attributes of events including case identifier,
activity name, and timestamp. Accordingly, we define a simple event, trace,
and event log. In the following, we introduce some basic concepts and
notations.

Let $A$ be a set. $A^{*}$ is the set of all finite sequences over $A$, and
$\mathcal{B}(A)$ is the set of all multisets over the set $A$. For
$A_{1},A_{2}\in\mathcal{B}(A)$, $A_{1}\subseteq A_{2}$ if for all $a\in A$,
$A_{1}(a)\leq A_{2}(a)$. A finite sequence over $A$ of length $n$ is a mapping
$\sigma\in\\{1,...,n\\}\rightarrow{A}$, represented as $\sigma=\langle
a_{1},a_{2},...,a_{n}\rangle$ where $\sigma_{i}=a_{i}=\sigma(i)$ for any
$1\leq i\leq n$, and $|\sigma|=n$. $a\in\sigma\Leftrightarrow{a=a_{i}}$ for
$1\leq i\leq n$. For $\sigma_{1},\sigma_{2}\in A^{*}$,
$\sigma_{1}\sqsubseteq\sigma_{2}$ if $\sigma_{1}$ is a subsequence of
$\sigma_{2}$, e.g., $\langle a,b,c,x\rangle\sqsubseteq\langle
z,x,a,b,b,c,a,b,c,x\rangle$. For $\sigma\in A^{*}$, $\\{a\in\sigma\\}$ is the
set of elements in $\sigma$, and $[a\in\sigma]$ is the multiset of elements in
$\sigma$, e.g., $[a\in\langle x,y,z,x,y\rangle]=[x^{2},y^{2},z]$.

######  Definition 1 (Simple Event)

A simple event is a tuple $e=(c,a,t)$, where $c\in\mathcal{C}$ is the case
identifier, $a\in\mathcal{A}$ is the activity associated to event $e$, and
$t\in\mathcal{T}$ is the timestamp of event $e$. $\pi_{X}(e)$ is the
projection of event $e$ on the attribute from domain $X$, e.g.,
$\pi_{\mathcal{A}}(e)=a$. We call
$\xi=\mathcal{C}\times\mathcal{A}\times\mathcal{T}$ the event universe.

######  Definition 2 (Simple Trace)

Let $\xi$ be the universe of events. A trace $\sigma=\langle
e_{1},e_{2},...,e_{n}\rangle$ in an event log is a sequence of events, i.e.,
$\sigma\in\xi^{*}$, s.t., for each $e_{i},e_{j}\in\sigma$:
$\pi_{\mathcal{C}}(e_{i})=\pi_{\mathcal{C}}(e_{j})$, and
$\pi_{\mathcal{T}}(e_{i})\leq\pi_{\mathcal{T}}(e_{j})$ if $i<j$. A simple
trace is a trace where all the events are projected on the activity attribute,
i.e., $\sigma\in\mathcal{A}^{*}$.

######  Definition 3 (Simple Event Log)

A simple event log is a multiset of simple traces, i.e.,
$L\in\mathcal{B}(\mathcal{A}^{*})$. We assume each trace in an event log
belongs to an individual and $\sigma\neq\langle\rangle$ if $\sigma\in L$.
$A_{L}=\\{a\in\mathcal{A}\mid\exists_{\sigma\in L}a\in\sigma\\}$ is the set of
activities in the event log $L$. $\tilde{L}=\\{\sigma\in L\\}$ is the set of
unique traces (variants) in the event log $L$. We denote $\mathcal{U}_{L}$ as
the universe of event logs.

######  Definition 4 (Trace Frequency)

Let $L$ be an event log, $f_{L}\in\tilde{L}\rightarrow{[0,1]}$ is a function
which retrieves the relative frequency of a trace in the event log $L$, i.e.,
$f_{L}(\sigma)=\nicefrac{{L(\sigma)}}{{|L|}}$ and
$\sum_{\sigma\in\tilde{L}}f_{L}(\sigma)=1$.

######  Definition 5 (Event Log Entropy)

$ent\in\mathcal{U}_{L}\rightarrow{\mathbb{R}_{\geq 0}}$ is a function which
retrieves the entropy of traces in an event log, s.t., for
$L\in\mathcal{U}_{L}$,
$ent(L)=-\sum_{\sigma\in\tilde{L}}f_{L}(\sigma)log_{2}{f_{L}(\sigma)}$. We
denote $max\\_ent(L)$ as the maximal entropy achieved when all the traces in
the event log are unique, i.e., $|\tilde{L}|=|L|$.

##  4 Privacy Quantification

We employ a risk-utility model for quantifying privacy in process mining where
disclosure risk and utility loss are measured to assess the effectiveness of
privacy preservation techniques before and after applying the techniques.

###  4.1 Disclosure Risk

In this subsection, we introduce identity/case disclosure and attribute/trace
disclosure for quantifying disclosure risk of event logs. Identity disclosure
quantifies how uniquely the trace owners, i.e., cases, can be re-identified.
Attribute disclosure quantifies how confidently the sensitive attributes of
cases (as individuals) can be specified. As discussed in [17], traces play the
role of both quasi-identifiers and sensitive attributes. That is, a complete
sequence of activities, which belongs to a case, is sensitive person-specific
information. At the same time, knowing a part of this sequence, as background
knowledge, can be exploited to re-identify the trace owner. In a simple event
log, traces, i.e., sequence of activities, are the only available information.
Therefore, attribute disclosure can be seen as trace disclosure.

In the following, we define set, multiset, and sequence as three types of
background knowledge based on traces in simple event logs that can be
exploited for uniquely re-identifying the trace owners or certainly specifying
their complete sequence of activities. Moreover, we consider a size for
different types of background knowledge as their power, e.g, the set
background knowledge of size 3 is more powerful than the same type of
background knowledge of size 2. Note that the assumed types of background
knowledge are the most general ones, and more types can be explored. However,
the general approach will be the same.

######  Definition 6 (Background Knowledge 1 - Set)

In this scenario, we assume that an adversary knows a subset of activities
performed for the case, and this information can lead to the identity or
attribute disclosure. Let $L$ be an event log, and $A_{L}$ be the set of
activities in the event log $L$. We formalize this background knowledge by a
function $proj_{set}^{L}\in 2^{A_{L}}\rightarrow{2^{L}}$. For $A\subseteq
A_{L}$, $proj_{set}^{L}(A)=[\sigma\in L\mid A\subseteq\\{a\in\sigma\\}]$. We
denote $cand_{set}^{l}(L)=\\{A\subseteq A_{L}\mid|A|=l\wedge
proj^{L}_{set}(A)\neq[]\\}$ as the set of all subsets over the set $A_{L}$ of
size $l$ for which there exists matching traces in the event log.

######  Definition 7 (Background Knowledge 2 - Multiset)

In this scenario, we assume that an adversary knows a sub-multiset of
activities performed for the case, and this information can lead to the
identity or attribute disclosure. Let $L$ be an event log, and $A_{L}$ be the
set of activities in the event log $L$. We formalize this background knowledge
by a function $proj_{mult}^{L}\in\mathcal{B}(A_{L})\rightarrow{2^{L}}$. For
$A\in\mathcal{B}(A_{L})$, $proj_{mult}^{L}(A)=[\sigma\in L\mid
A\subseteq[a\in\sigma]]$. We denote
$cand_{mult}^{l}(L)=\\{A\in\mathcal{B}(A_{L})\mid|A|=l\wedge
proj^{L}_{mult}(A)\neq[]\\}$ as the set of all sub-multisets over the set
$A_{L}$ of size $l$ for which there exists matching traces in the event log.

######  Definition 8 (Background Knowledge 3 - Sequence)

In this scenario, we assume that an adversary knows a subsequence of
activities performed for the case, and this information can lead to the
identity or attribute disclosure. Let $L$ be an event log, and $A_{L}$ be the
set of activities in the event log $L$. We formalize this background knowledge
by a function $proj_{seq}^{L}\in A_{L}^{*}\rightarrow{2^{L}}$. For $\sigma\in
A_{L}^{*}$, $proj_{seq}^{L}(\sigma)=[\sigma^{\prime}\in
L\mid\sigma\sqsubseteq\sigma^{\prime}]$. We denote
$cand_{seq}^{l}(L)=\\{\sigma\in A_{L}^{*}\mid|\sigma|=l\wedge
proj^{L}_{seq}(\sigma)\neq[]\\}$ as the set of all subsequences of size
(length) $l$, based on the activities in $A_{L}$, for which there exists
matching traces in the event log.

######  Example 1 (background knowledge)

Let $L=[\langle a,b,c,d\rangle^{10},\langle a,c,b,d\rangle^{20},\langle
a,d,\\\ b,d\rangle^{5},\langle a,b,d,d\rangle^{15}]$ be an event log.
$A_{L}=\\{a,b,c,d\\}$ is the set of unique activities, and
$cand_{set}^{2}(L)=\\{\\{a,b\\},\\{a,c\\},\\{a,d\\},\\{b,c\\},\\{b,d\\},\\{d,c\\}\\}$
is the set of candidates of the set background knowledge of size $2$. For
$A=\\{b,d\\}\in cand_{set}^{2}(L)$ as a candidate of the set background
knowledge of size $2$, $proj_{set}^{L}(A)=[\langle a,b,c,d\rangle^{10},\langle
a,c,b,d\rangle^{20},\langle a,d,b,d\rangle^{5},\langle a,b,d,d\rangle^{15}]$.
For $A=[b,d^{2}]$ as a candidate of the multiset background knowledge,
$proj_{mult}^{L}(A)=[\langle a,d,b,d\rangle^{5},\langle a,b,d,d\rangle^{15}]$.
Also, for $\sigma=\langle b,d,d\rangle$ as a candidate of the sequence
background knowledge, $proj_{seq}^{L}(\sigma)=[\langle a,b,d,d\rangle^{15}]$.

As Example 1 shows, the strength of background knowledge from the weakest to
the strongest w.r.t. the type is as follows: set, multiset, and sequence,
i.e., given the event log $L$, $proj^{L}_{seq}(\langle b,d,d\rangle)\subseteq
proj^{L}_{mult}([b,d^{2}])\subseteq proj^{L}_{set}(\\{b,d\\})$.

####  4.1.1 Identity (Case) Disclosure

We use the uniqueness of traces w.r.t. the background knowledge of size $l$ to
measure the corresponding case disclosure risk in an event log. Let $L$ be an
event log and $type\in\\{set,mult,seq\\}$ be the type of background knowledge.
The case disclosure based on the background knowledge $type$ of size $l$ is
calculated as follows:

$\small cd_{type}^{l}(L)=\sum_{x\in
cand_{type}^{l}(L)}\frac{\nicefrac{{1}}{{|proj^{L}_{type}(x)|}}}{|cand_{type}^{l}(L)|}$
(1)

Equation (1) calculates the average uniqueness based on the candidates of
background knowledge, i.e., $x\in cand_{type}^{l}(L)$. Note that we consider
equal weights for the candidates of background knowledge. However, they can be
weighted based on the various criteria, e.g., the sensitivity of the
activities included. One can also consider the worst case, i.e., the maximal
uniqueness, rather than the average value.

######  Example 2 (insufficiency of case disclosure analysis)

Consider $L_{1}=[\langle a,\\\ b,c,d\rangle,\langle a,c,b,d\rangle,\langle
a,b,c,c,d\rangle,\langle a,b,b,c,d\rangle]$ and $L_{2}=[\langle
a,b,c,d\rangle^{4},\langle e,f\rangle^{4},\langle g,h\rangle^{4}]$ as two
event event logs. $A_{L_{1}}=\\{a,b,c,d\\}$ and
$A_{L_{2}}=\\{a,b,c,d,e,f,g,h\\}$ are the set of unique activities in $L_{1}$
and $L_{2}$, respectively.
$cand_{set}^{1}(L_{1})=\\{\\{a\\},\\{b\\},\\{c\\},\\{d\\\ \\}\\}$ and
$cand_{set}^{1}(L_{2})=\\{\\{a\\},\\{b\\},\\{c\\},\\{d\\},\\{e\\},\\{f\\},\\{g\\},\\{h\\}\\}$
are the set of candidates of the set background knowledge of size $1$. Both
event logs have the same value as the case disclosure for the set background
knowledge of size $1$ ($cd_{set}^{1}(L_{1})=cd_{set}^{1}(L_{2})=1/4$).
However, in $L_{2}$, the complete sequence of activities performed for a
victim case is disclosed by knowing only one activity without uniquely
identifying the corresponding trace.

Example 2 clearly shows that measuring the uniqueness alone is insufficient to
demonstrate disclosure risks in event logs and the uncertainty in the set of
sensitive attributes matching with the assumed background knowledge need to be
measured, as well. In the following, we define a measure to quantify the
uncertainty in the set of matching traces. Note that, the same approach can be
exploited to quantify the disclosure risk of any other sensitive attribute
matching with some background knowledge.

####  4.1.2 Attribute (Trace) Disclosure

We use the entropy of matching traces w.r.t. background knowledge of size $l$
to measure the corresponding trace disclosure risk in an event log. Let $L$ be
an event log and $type\in\\{set,mult,seq\\}$ be the type of background
knowledge. The trace disclosure based on the background knowledge $type$ of
size $l$ is calculated as follows:

$\small td_{type}^{l}(L){=}1-{\sum_{x\in
cand_{type}^{l}(L)}{\frac{\nicefrac{{ent(proj^{L}_{type}(x))}}{{max\\_ent(proj^{L}_{type}(x))}}}{|cand_{type}^{l}(L)|}}}$
(2)

In (2), $max\\_ent(proj^{L}_{type}(x))$ is the maximal entropy for the
matching traces based on the type and size of background knowledge, i.e.,
uniform distribution of the matching traces. As discussed for (1), in (2), we
also assume equal weights for the candidates of background knowledge. However,
one can consider different weights for the candidates. Also, the worst case,
i.e., the minimal entropy, rather than the average entropy can be considered.

The trace disclosure of the event logs in Example 2 is as follows:
$td_{set}^{1}(L_{1})=0$ (the multiset of matching traces has the maximal
entropy) and $td_{set}^{1}(L_{2})=1$ (the entropy of matching traces is 0).
These results distinguish the disclosure risk of the event logs.

###  4.2 Utility Loss

In this subsection, we introduce a measure based on the earth mover’s distance
[19] for quantifying the utility loss after applying a privacy preservation
technique to an event log. The earth mover’s distance describes the distance
between two distributions. In an analogy, given two piles of earth, it
expresses the effort required to transform one pile into the other. First, we
introduce the concept of reallocation indicating how an event log is
transformed into another event log. Then, we define a trace distance function
expressing the cost of transforming one trace into another one. Finally, we
introduce the utility loss measure that indicates the entire cost of
transforming an event log to another one using the introduced reallocation and
distance functions.

####  4.2.1 Reallocation

Let $L$ be the original event log and $L^{\prime}$ be an anonymized event log
derived from the original event log. We introduce
$r\in\tilde{L}\times\tilde{L}^{\prime}\rightarrow{[0,1]}$ as a function that
indicates the movement of frequency between two event logs.
$r(\sigma,\sigma^{\prime})$ describes the relative frequency of
$\sigma\in\tilde{L}$ that should be transformed to
$\sigma^{\prime}\in\tilde{L}^{\prime}$. To make sure that a reallocation
function properly transforms $L$ into $L^{\prime}$, the frequency of each
$\sigma\in\tilde{L}$ should be considered, i.e., for all $\sigma\in\tilde{L}$,
$f_{L}(\sigma)=\sum_{\sigma^{\prime}\in\tilde{L}^{\prime}}r(\sigma,\sigma^{\prime})$.
Similarly, the probability mass of traces
$\sigma^{\prime}\in\tilde{L}^{\prime}$ should be preserved, i.e., for all
$\sigma^{\prime}\in\tilde{L}^{\prime}$,
$f_{L^{\prime}}(\sigma^{\prime})=\sum_{\sigma\in\tilde{L}}r(\sigma,\sigma^{\prime})$.
We denote $\mathcal{R}$ as the set of all reallocation functions which depends
on $L$ and $L^{\prime}$.

Table 1: The dissimilarity between two event logs based on the earth mover’s distance assuming $r_{s}$ as a reallocation function and $d_{s}$ as the normalized Levenshtein distance. $r_{s}\cdot d_{s}$ | $\langle a,b,c,d\rangle$ | $\langle a,c,b,d\rangle$ | $\langle a,e,c,d\rangle^{49}$ | $\langle a,e,b,d\rangle^{49}$
---|---|---|---|---
$\langle a,b,c,d\rangle^{50}$ | $0.01\cdot 0$ | $0\cdot 0.5$ | $0.49\cdot 0.25$ | $0\cdot 0.5$
$\langle a,c,b,d\rangle^{50}$ | $0\cdot 0.5$ | $0.01\cdot 0$ | $0\cdot 0.5$ | $0.49\cdot 0.25$

####  4.2.2 Trace Distance

A trace distance function
$d\in\mathcal{A}^{*}\times\mathcal{A}^{*}\rightarrow{[0,1]}$ expresses the
distance between traces. This function is 0 if and only if two traces are
equal, i.e., $d(\sigma,\sigma^{\prime})=0{\iff}\sigma=\sigma^{\prime}$. This
function should also be symmetrical, i.e.,
$d(\sigma,\sigma^{\prime})=d(\sigma^{\prime},\sigma)$. Different distance
functions can be considered satisfying these conditions. We use the normalized
string edit distance (Levenshtein) [9].

####  4.2.3 Utility Loss

Let $L$ be an original event log, and $L^{\prime}$ be an anonymized event log
derived from the original event log. Several reallocation functions might
exist. However, the earth mover’s distance problem aims to express the
shortest distance between the two event logs, i.e., the least mass movement
over the least distance between traces. Therefore, the difference between $L$
and $L^{\prime}$ using a reallocation function $r$ is the inner product of
reallocation and distance. The data utility preservation is considered as
$du(L,L^{\prime})=1-\min\limits_{r\in\mathcal{R}}ul(r,L,L^{\prime})$.

$\small ul(r,L,L^{\prime})=r\cdot
d=\sum_{\sigma\in\tilde{L}}\sum_{\sigma^{\prime}\in\tilde{L}^{\prime}}{r(\sigma,\sigma^{\prime})d(\sigma,\sigma^{\prime})}$
(3)

###### Example 3

(using earth mover’s distance to calculate dissimilarity between event logs)
Let $L=[\langle a,b,c,d\rangle,\langle a,c,b,d\rangle,\langle
a,e,c,d\rangle^{49},\langle a,e,b,d\rangle^{49}]$ and $L^{\prime}=[\langle
a,b,c,d\rangle^{50},\langle a,c,b,d\rangle^{50}]$ be the original and
aninymized event logs, respectively. Table 1 shows the calculations assuming
$r_{s}$ as a reallocation function and $d_{s}$ as the normalized Levenshtein
distance, e.g., $r_{s}(\langle a,b,c,d\rangle,\langle a,e,c,d\rangle)=0.49$
and $d_{s}(\langle a,b,c,d\rangle,\langle a,e,c,d\rangle)=0.25$.
$ul(r_{s},L,L^{\prime})=0.24$ and $du(L,L^{\prime})=0.76$.

Table 2: The general statistics of the event logs used in the experiments. Event Log | #traces | #variants | #events |  | #unique_activities
---
| trace_uniqueness
---
Sepsis-Cases [10] | 1050 | 845 | 15214 | 16 | 80%
BPIC-2017-APP [21] | 31509 | 102 | 239595 | 10 | 0.3%

##  5 Experiments

In this section, we demonstrate the experiments on real-life event logs to
advocate the proposed measures. We employ two human-centered event logs, where
the case identifiers refer to individuals. Sepsis-Cases [10] is a real-life
event log containing events of sepsis cases from a hospital. BPIC-2017-APP
[21] is also a real-life event log pertaining to a loan application process of
a Dutch financial institute. We choose these event logs because they are
totally different w.r.t. the uniqueness of traces. Table 2 shows the general
statistics of these event logs. Note that variants are the unique traces, and
$trace\\_uniquness=\nicefrac{{\\#variants}}{{\\#traces}}$. The implementation
as a Python program is available on
GitHub.111https://github.com/m4jidRafiei/privacy_quantification

###  5.1 Disclosure Risk Analysis

In this subsection, we show the functionality of the proposed measures for
disclosure risk analysis. To this end, we consider three types of background
knowledge (set, multiset, and sequence) and vary the background knowledge
power (size) from 1 to 6. Figure 2a shows the results for the Sepsis-Cases
event log where the uniqueness of traces is high. As shown, the disclosure
risks are higher for the more powerful background knowledge w.r.t. the type
and size.

Figure 2b demonstrates the results for the BPIC-2017-APP event log, where the
uniqueness of traces is low. As shown, the case disclosure risk is low, which
is expected regarding the low uniqueness of traces. However, the trace
disclosure risk is high which indicates low entropy (uncertainty) of the
traces. Moreover, for the stronger background knowledge w.r.t. the size, one
can assume a higher case disclosure risk. However, the trace disclosure risk
is correlated with the entropy of the sensitive attribute values and can be a
high value even for weak background knowledge. The above-mentioned analyses
clearly show that uniqueness alone cannot reflect the actual disclosure risk
in an event log.

(a) Sepsis-Cases [10].

(b) BPIC-2017-APP [21].

Figure 2: Analyses of the case disclosure ($cd$) and the trace disclosure
($td$) based on the three types of background knowledge (i.e., $set$, $mult$,
and $seq$) when we vary the background knowledge power (size) from 1 to 6. For
example, in the Sepsis-Cases event log, the case disclosure risk of the
background knowledge $seq$ ($cd\\_seq$) of size 3 is 0.188.

###  5.2 Utility Loss Analysis

In this subsection, we demonstrate the functionality of the proposed measure
in Section 4.2 for quantifying data utility preservation after applying a
privacy preservation technique. We use PPDP-PM [15] as a privacy preservation
tool for process mining to apply the $\mathrm{TLKC}$-privacy model [17] to a
given event log. The $\mathrm{TLKC}$-privacy model is a group-based privacy
preservation technique which provides a good level of flexibility through
various parameters such as the type and size (power) of background knowledge.
The $\mathrm{T}$ in this model refers to the accuracy of timestamps in the
privacy-aware event log, $\mathrm{L}$ refers to the power of background
knowledge222Note that this $\mathrm{L}$ is identical to the $l$ introduced as
the power (size) of background knowledge and should not be confused with $L$
as the event log notation., $\mathrm{K}$ refers to the $k$ in the
$k$-anonymity definition [20], and $\mathrm{C}$ refers to the bound of
confidence regarding the sensitive attribute values in an equivalence class.

(a) Using set as background knowledge.

(b) Using sequence as background knowledge.

Figure 3: The utility loss and disclosure risk analyses for the Sepsis-Cases event log where the background knowledge is set or sequence, and the power (size) of background knowledge is 6. Table 3: The general statistics before and after applying the $\mathrm{TLKC}$-privacy model. Event Log | #traces | #variants | #events | #unique_activities
---|---|---|---|---
Original Sepsis-Cases | 1050 | 845 | 15214 | 16
Anonymized Sepsis-Cases |  BK type=set BK size (L)=6 | K=20 | 1050 | 842 | 15103 | 12
K=40 | 1050 | 842 | 14986 | 11
K=60 | 1050 | 818 | 14809 | 11
BK type=seq BK size (L)=6 | K=20 | 1050 | 34 | 3997 | 6
K=40 | 1050 | 33 | 4460 | 5
K=60 | 1050 | 18 | 3448 | 4

Assuming set (6) and sequence (8) as the types of background knowledge, we
apply the $\mathrm{TLKC}$-privacy model to the Sepsis-Cases event log with the
following parameters: $\mathrm{L}=6$ (as the maximum background knowledge
power in our experiments), $\mathrm{K}=\\{20,40,60\\}$, $\mathrm{C}=1$ (there
is no additional sensitive attribute in a simple event log), and $\mathrm{T}$
is set to the maximal precision ($\mathrm{T}$ has no effect on a simple event
log). That is, the $\mathrm{TLKC}$-privacy model is simplified to
$k$-anonymity where the quasi-identifier (background knowledge) is the set or
sequence of activities. Table 3 demonstrates the general statistics of the
event logs before and after applying the privacy preservation technique.

Figure 3a shows disclosure risk and data utility analyses for the background
knowledge set, and Fig. 3b shows the same analyses for the background
knowledge sequence. In both types of background knowledge, the data utility
value decreases. For the stronger background knowledge, i.e., sequence, the
utility loss is much higher which is expected w.r.t. the general statistics in
Table 3. However, the data utility for the weaker background knowledge remains
high which again complies with the general statistics. Note that since we
apply $k$-anonymity (simplified $\mathrm{TLKC}$-model) only case disclosure,
which is based on the uniqueness of traces, decreases. Moreover, for the
sequence background knowledge, higher values for $\mathrm{K}$ result in more
similar traces. Therefore, the trace disclosure risk, in the anonymized event
logs, drastically increases. These analyses demonstrate that privacy
preservation techniques should consider different aspects of disclosure risk
while balancing data utility preservation and sensitive data protection.

##  6 Conclusion

Event logs often contain highly sensitive information, and regarding the rules
imposed by regulations, these sensitive data should be analyzed responsibly.
Therefore, privacy preservation in process mining is recently receiving more
attention. Consequently, new measures need to be defined to evaluate the
effectiveness of the privacy preservation techniques both from the sensitive
data protection and data utility preservation point of views. In this paper,
using a trade-off approach, we introduced two measures for quantifying
disclosure risks: identity/case disclosure and attribute/trace disclosure, and
one measure for quantifying utility loss. The introduced measures were applied
to two real-life event logs. We showed that even simple event logs could
reveal sensitive information. Moreover, for the first time, the effect of
applying a privacy preservation technique on data utility rather than result
utility was explored. The data utility measure is based on the earth mover’s
distance and can be extended to evaluate the utility w.r.t. the different
perspectives of process mining, e.g., time, resource, etc.

## Acknowledgment

Funded under the Excellence Strategy of the Federal Government and the Länder.
We also thank the Alexander von Humboldt (AvH) Stiftung for supporting our
research.

## References

  * [1] van der Aalst, W.M.P.: Process Mining - Data Science in Action, Second Edition. Springer (2016). https://doi.org/10.1007/978-3-662-49851-4
  * [2] van der Aalst, W.M.P.: Responsible data science: using event data in a “people friendly” manner. In: International Conference on Enterprise Information Systems. pp. 3–28. Springer (2016)
  * [3] Bertino, E., Fovino, I.N., Provenza, L.P.: A framework for evaluating privacy preserving data mining algorithms. Data Min. Knowl. Discov. 11(2), 121–154 (2005)
  * [4] Bertino, E., Lin, D., Jiang, W.: A survey of quantification of privacy preserving data mining algorithms. In: Aggarwal, C.C., Yu, P.S. (eds.) Privacy-Preserving Data Mining - Models and Algorithms, Advances in Database Systems, vol. 34, pp. 183–205. Springer (2008)
  * [5] Elkoumy, G., Fahrenkrog-Petersen, S.A., Dumas, M., Laud, P., Pankova, A., Weidlich, M.: Secure multi-party computation for inter-organizational process mining. In: Enterprise, Business-Process and Information Systems Modeling - 21st International Conference, BPMDS. Springer (2020)
  * [6] Fahrenkrog-Petersen, S.A., van der Aa, H., Weidlich, M.: PRETSA: event log sanitization for privacy-aware process discovery. In: International Conference on Process Mining, ICPM 2019, Aachen, Germany (2019)
  * [7] Iyengar, V.S.: Transforming data to satisfy privacy constraints. In: Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 279–288. ACM (2002)
  * [8] Jr., R.J.B., Agrawal, R.: Data privacy through optimal k-anonymization. In: Proceedings of the 21st International Conference on Data Engineering, ICDE (2005)
  * [9] Levenshtein, V.I.: Binary codes capable of correcting deletions, insertions, and reversals. In: Soviet physics doklady. vol. 10, pp. 707–710 (1966)
  * [10] Mannhardt, F.: Sepsis cases-event log. Eindhoven University of Technology (2016)
  * [11] Mannhardt, F., Koschmider, A., Baracaldo, N., Weidlich, M., Michael, J.: Privacy-preserving process mining - differential privacy for event logs. Business & Information Systems Engineering 61(5), 595–614 (2019)
  * [12] Michael, J., Koschmider, A., Mannhardt, F., Baracaldo, N., Rumpe, B.: User-centered and privacy-driven process mining system design for IoT. In: Information Systems Engineering in Responsible Information Systems. pp. 194–206 (2019)
  * [13] Pika, A., Wynn, M.T., Budiono, S., ter Hofstede, A.H., van der Aalst, W.M.P., Reijers, H.A.: Privacy-preserving process mining in healthcare. International Journal of Environmental Research and Public Health 17(5), 1612 (2020)
  * [14] Rafiei, M., van der Aalst, W.M.P.: Mining roles from event logs while preserving privacy. In: Business Process Management Workshops - BPM 2019 International Workshops, Vienna, Austria. pp. 676–689 (2019)
  * [15] Rafiei, M., van der Aalst, W.M.P.: Practical aspect of privacy-preserving data publishing in process mining. In: Proceedings of the Best Dissertation Award, Doctoral Consortium, and Demonstration & Resources Track at BPM 2020 co-located with the 18th International Conference on Business Process Management (BPM 2020). CEUR-WS.org (2020)
  * [16] Rafiei, M., van der Aalst, W.M.P.: Privacy-preserving data publishing in process mining. In: Business Process Management Forum - BPM Forum 2020, Seville, Spain, September 13-18. pp. 122–138. Springer (2020)
  * [17] Rafiei, M., Wagner, M., van der Aalst, W.M.P.: TLKC-privacy model for process mining. In: Research Challenges in Information Science - 14th International Conference, RCIS. pp. 398–416. Springer International Publishing (2020)
  * [18] Rafiei, M., von Waldthausen, L., van der Aalst, W.M.P.: Supporting condentiality in process mining using abstraction and encryption. In: Data-Driven Process Discovery and Analysis - 8th IFIP WG 2.6 International Symposium, SIMPDA 2018, and 9th International Symposium, SIMPDA 2019, Revised Selected Papers (2019)
  * [19] Rüschendorf, L.: The wasserstein distance and approximation theorems. Probability Theory and Related Fields 70(1), 117–129 (1985)
  * [20] Sweeney, L.: k-anonymity: A model for protecting privacy. International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems 10(05), 557–570 (2002)
  * [21] Van Dongen, B.F.: BPIC 2017. Eindhoven University of Technology (2017)
  * [22] von Voigt, S.N., Fahrenkrog-Petersen, S.A., Janssen, D., Koschmider, A., Tschorsch, F., Mannhardt, F., Landsiedel, O., Weidlich, M.: Quantifying the re-identification risk of event logs for process mining - empiricial evaluation paper. In: Advanced Information Systems Engineering, CAiSE (2020)