# Using stochastic resonance of two-level systems to increase qubit
decoherence times

Yujun Choi  Department of Physics, Virginia Tech, Blacksburg, VA 24061, United
States  Virginia Tech Center for Quantum Information Science and Engineering,
Blacksburg, VA 24061, United States  S. N. Coppersmith  School of Physics, The
University of New South Wales, Sydney, NSW 2052, Australia  Robert Joynt
<EMAIL_ADDRESS>Department of Physics, University of Wisconsin-Madison,
Madison, WI 53706, United States  Theoretical Sciences Visiting Program
(TSVP), Okinawa Institute of Science and Technology Graduate University, Onna,
904-0495, Japan

###### Abstract

Two-level systems (TLS) are the major source of dephasing of spin qubits in
numerous quantum computing platforms. In spite of much effort, it has been
difficult to substantially mitigate the effects of this noise or, in many
cases, to fully understand its physical origin. We propose a method to make
progress on both of these issues. When an oscillating field is applied to a
TLS, stochastic resonance can occur and the noise spectrum is moved to higher
frequencies. This shift in the TLS noise spectrum will increase the dephasing
times of the qubits that they influence. Furthermore, the details of this
effect depend on the physical properties of the noise sources. Thus one can
use qubit spectroscopy to investigate their physical properties, specifically
the extent to which the TLS themselves possess quantum coherence. We find that
it should be possible to determine the dephasing rate and the energy level
separation of the TLS themselves in this way.

##  I Introduction

$1/f$ noise is a major source of decoherence in many quantum computing
platforms [1]. It is generally believed to originate from two-level systems
(TLS) whose fluctuations produce a random time-dependent field at the qubit
[2, 3, 4, 5, 6, 7]. To reduce the effect of $1/f$ noise, several experimental
techniques have been proposed. By far the best known is dynamical decoupling
(DD) that compensates accumulated phase error by applying refocusing pulses to
a qubit [8, 9]. DD can substantially increase the dephasing time of the qubit.
However, it has the disadvantages of introducing more hardware and additional
manipulation of the qubit.

Instead of driving the qubits to increase decoherence times, one may instead
drive the noise sources that affect them. We shall call this noise source
driving (NSD). It is much less studied than dynamical decoupling. One example
of NSD is spin bath driving which has shown to be useful for nitrogen-vacancy
center spin qubits in diamond. This scheme applies a radio-frequency pulse or
continuous wave to surface spins on the diamond, giving rise to partial
decoupling of the dipolar interaction to the NV spin qubits [10, 11]. The spin
bath driving shifts the noise power spectral density (PSD) from low to high
frequencies, which increases the dephasing times [12, 13]. Experiment shows
that this works well with both monochromatic and polychromatic driving fields
[14]. This method is somewhat analogous to dynamical decoupling but it has the
advantage that one does not manipulate the qubits - the carriers of quantum
information - but rather the noise sources, which carry no useful information
for quantum computing. Hence we would expect that the technical requirements
on the source of the driving fields are much less stringent for NSD than for
DD.

Another NSD scheme for superconducting qubits has also been proposed.
Theoretical study shows that the coherence times of a superconducting qubit
can be increased by adding noise to the surrounding TLSs [15]. It is
noteworthy that the relaxation (decoherence) time of the qubit is closely
related to the decoherence (relaxation) of the TLSs. An experimental study
exhibits a 23% increase in the relaxation time of a transmon qubit on average
when a DC electric field is applied to the TLSs [16].

In the present work, we attempt to generalize and systematize the theory of
NSD by using concepts from the established field of stochastic resonance (SR).
SR as originally conceived is a phenomenon in which a charged particle moving
back and forth randomly in a bistable potential is driven by an applied
periodic oscillating electromagnetic field [17]. It is certainly also possible
to drive the system with sound waves, and this may be more practical in some
circumstances [18, 15]. In the absence of driving, the dynamics are purely
relaxational, typically leading to a PSD that is Lorentzian. When driven, one
sees a modification of the PSD of the particle’s motion that depends on the
relation of the period of the oscillating field with the intrinsic relaxation
rate. SR differs from ordinary resonance in that it does not require that the
system in question has a restoring force. Nevertheless, the modifications
produced in the PSD can be dramatic, and this has important consequences for
decoherence.

The theory of SR was developed originally in the weak-field limit [19] and
later in the strong-field limit [20, 21]. A more general approach that is
independent of external signal strength was also studied [22]. In recent years
quantum stochastic resonance was observed experimentally in the tunneling
process of electrons in a quantum dot [23] and in the spin dynamics of a
single Fe atom under scanning tunneling microscope [24]. The movement of
weight to higher frequencies is maximized when the driving frequency is about
half of the switching rate of the TLS [17].

All of the above-cited theoretical work assumed essentially classical dynamics
for the TLS. The extension to quantum stochastic resonance where a spin-boson
model was adopted instead of a classical particle in a bistable potential has
also been examined [25, 26]. The traditional assumption is that relaxation
times from the phonon bath are so short compared to the tunneling rate that it
is justified to adopt a Markovian picture in which the transition rates do not
reflect quantum coherence and therefore have no significant time dependence.
However, the actual evidence for this is slim, and of course, it may be true
in some systems and not in others. Hence it is important to have signatures of
quantumness that can only be obtained using a theoretical framework in which
coherent quantum effects can be present but which yet has a well-defined
classical limit. The Lindblad equation satisfies these criteria, which is why
we use it in this paper.

This paper has two goals. The first is to show that the shift in the PSD that
SR produces can mitigate qubit decoherence, even in the presence of quantum
coherence of the TLS. The second is to show that NSD can reveal how much
quantum coherence a TLS has, and of what type.

In Sec. II we review the theoretical background for SR. In Sec. III, the model
for quantum TLS is introduced. We discuss the meaning of the parameters in the
model and give some computational details. In Sec. IV we compute the shift in
the noise PSD for a single TLS due to the driving and its consequences for the
qubit dephasing time $T_{\phi}$. In Sec. V we deal with the case where the TLS
dynamics are not purely relaxational. In Sec. VI we extend this to multiple
noise sources. Finally, in Sec. VII we summarize the results for the
mitigation of qubit decoherence and discuss the question of how much light NSD
experiments can throw on the quantum nature of the noise sources.

##  II Stochastic resonance and Qubit Dephasing

Classical SR is a phenomenon that occurs when a classical two-state system is
driven by a sinusoidally oscillating external field. The dynamics of the
undriven system are assumed to be purely relaxational so its PSD has a single
peak at zero frequency. When driven, the system also oscillates at higher
frequencies and the height of the central peak decreases.

###  II.1 Classical Model

We establish notation and conceptual background in this section. The treatment
follows Refs. [19] and [22]. The two states of the system, (henceforth just
the “TLS”) are denoted by $s(t)=\pm 1$. If at time $t$ we have that $s(t)=\pm
1$, then at time $t+dt$ the system has made a transition to $s(t)=\mp 1$ with
probability $W_{\pm}(t)dt$. We denote the probability that the TLS is in the
state $\pm 1$ state at time $t$ by $P(\pm 1,t)$. Since $P(+1,t)+P(-1,t)=1$ it
makes sense to define $s_{z}(t)=P(+1,t)-P(-1,t)$. The differential equation
for $s_{z}(t)$ is

$\frac{ds_{z}(t)}{dt}=-W(t)s_{z}(t)+\delta W(t).$ (1)

Here $W=W_{-}+W_{+}$ and $\delta W=W_{-}-W_{+}$. This equation can be solved
exactly using the integrating factor method. Our main interest is in the
correlation function $S(\tau)=\langle s(t+\tau)s(t)\rangle$, where the angle
brackets denote an average over initial conditions and over $t$. The PSD is

$S(\omega)=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\tau
e^{-i\omega\tau}S(\tau).$ (2)

We have that $S(-\omega)=S(\omega)$ so we will usually plot $S(\omega)$ only
for $\omega>0$. The key point for the PSD is that since $S(\tau=0)=\langle
s_{z}^{2}(t)\rangle=1$ we find that

$\int_{-\infty}^{\infty}d\omega S(\omega)=S(\tau=0)=1.$ (3)

Thus the integral of the PSD is independent of the form of $W_{\pm}(t)$. This
implies that driving the TLS can move weight around in the integral of the PSD
but the total weight (the value of the integral) is unaffected.

If there is no driving then $W_{\pm}(t)=W_{\pm}$ is independent of time and we
find

$S(\omega)=\frac{2W/\pi}{W^{2}+\omega^{2}}$ (4)

in the simple case where $W_{+}=W_{-}=W/2$. We have set $\delta W=0$ here
since it only affects the average value of $s(t)$ which is less interesting
for our purposes.

To drive the TLS we choose a driving frequency $\omega_{d}>0$ and a driving
strength $A$ and then set $W_{\pm}=W/2+A\cos\omega_{d}t$. One then finds [19]

$S(\omega)=(1-R)\frac{2W/\pi}{W^{2}+\omega^{2}}+\frac{R}{2}[\delta(\omega-\omega_{d})+\delta(\omega+\omega_{d})].$
(5)

with

$R=\frac{A^{2}/2}{W^{2}+\omega_{d}^{2}}$ (6)

$R$ is the redistribution ratio that quantifies how much of the weight in the
PSD has been moved from low frequencies of order $W$ to frequencies greater
than $W$ (when $\omega_{d}>W$). The sum rule for $S(\omega)$ implies that
$0\leq R\leq 1$.

When choosing $\omega_{d}$ there is a trade-off. Large $\omega_{d}$ moves the
delta-function at $\omega=\omega_{d}$ further to the right in a plot of
$S(\omega)$ versus $\omega$, but the amplitudes of the delta functions
decrease as $[\omega_{d}^{2}+W^{2}]^{-1}$.

At higher order in $A$, harmonics appear in the PSD; there are additional
delta functions centered at $\omega_{d},2\omega_{d},3\omega_{d},...$, and the
noise spectrum is also modified [22]. In the Supplemental Material we show
that when $\omega_{d}\gg W$, the functional form of the noise spectrum does
not change significantly even when $A$ is large, so the $R$ ratio provides a
good characterization of the shift in noise power induced by the driving.

###  II.2 Effect on Qubit Dephasing

We now regard the TLS as a source of noise that acts on qubits. Our model of
the TLS is very general, and capable of describing magnetic noise, for
example. For concreteness, however, let us suppose that the TLS is a
fluctuating electric dipole $\mathbf{p}(t)$ so that the two states are
$\pm\mathbf{p}=\mathbf{p}_{0}s(t)$. Let the TLS be at the origin of
coordinates and a charge qubit be at $\mathbf{r}$. The operating frequency of
the qubit is proportional to $E_{i}(\mathbf{r})$, the $i$th component of the
electric field at the qubit. The field of the TLS at $\mathbf{r}$ is
$E_{i}(t)=-\partial V(t)/\partial x_{i}$. with
$V(t)=\mathbf{p}_{0}\cdot\mathbf{r}\,s(t)/4\pi\epsilon_{0}\epsilon_{r}r^{2}$,
in a medium with relative permittivity $\epsilon_{r}$. The dephasing time is
determined by the correlation function

$S_{E}(\tau)=\langle E_{i}(t+\tau)E_{i}(t)\rangle\propto\langle
s(t+\tau)s(t)\rangle$ (7)

which leads immediately to $S_{E}(\omega)\propto S(\omega)$. For some further
details about the application to a charge qubit, see the Supplemental
Material.

We shall, for the most part, focus in this paper on qubits whose dominant
decoherence channel is dephasing from low-frequency noise, and the associated
decoherence time will simply be called $T_{\phi}$. As long as $T_{\phi}\ll
T_{1}$, we also have $T_{\phi}=T_{2}$.

$S(\omega)$ controls dephasing as follows. The off-diagonal elements of the
qubit density matrix decay at short times according to
$\exp(-t^{2}/T_{\phi}^{2})$ with

$\frac{1}{T_{\phi}^{2}}=C_{0}\int_{0}^{\infty}S(\omega)F(\omega
T_{\phi})d\omega.$ (8)

Here $F$ is a window function that selects for the frequencies in $S(\omega)$
that satisfy $\omega T_{\phi}\leq 1$. For Gaussian noise we have $F(\omega
T_{\phi})=\sin^{2}(\omega T_{\phi}/2)/(\omega T_{\phi}/2)^{2}$. This is the
form of $F$ that we use below to compute PSDs. But for any kind of noise, $F$
has the same qualitative behavior. An approximation that captures the
essential physics is obtained by taking $F$ to be a step function, resulting
in

$\frac{1}{T_{\phi}^{2}}=C_{0}\int_{0}^{1/T_{\phi}}S(\omega)d\omega.$ (9)

$T_{\phi}$ is limited only by the low-frequency part of the noise. Of course,
this is a very crude way to compute $T_{\phi}$ but in the absence of more
details about the TLS there is little point in attempting to refine it. The
coupling constant $C_{0}$ depends on how the noise field is produced by the
TLS (as just described for the example of an electric dipole), and how
strongly the qubit couples to that field. Some additional details for a charge
qubit are again given in the Supplemental Material.

In Fig. 1 we give the graphical picture of how SR enhances $T_{\phi}$, using
9. In the absence of driving the noise is in the central Lorentzian and the
integral of that (together with $C_{0}$) gives $T_{\phi}$. When driven, the
weight in the PSD moves above $1/T_{\phi}$, which is the upper cutoff in the
integral, so $T_{\phi}$ is enhanced.

Figure 1: Solid line is the PSD in the absence of driving, while the dashed
line is the PSD of the driven TLS. The peak at $\omega_{d}$ represents a
$\delta$-function. When $1/T_{\phi}>W$ and the TLS is driven at a frequency
$\omega_{d}>1/T_{\phi}$, then the weight in the integral in Eq. 8 to the left
of the cutoff $1/T_{\phi}$ is reduced and $T_{\phi}$ increases.

When the inequalities $W\ll 1/T_{\phi}\ll\omega_{d}$ hold so that the
characteristic frequencies in the problem are well separated, then we can make
definite statements about the enhancement of $T_{\phi}$. Denoting the
dephasing time in the absence (presence) of driving as $T_{\phi,0}$
($T_{\phi,d}$), Eq. 9 leads to

$\frac{T_{\phi,d}}{T_{\phi,0}}=\frac{1}{\sqrt{1-R}}.$ (10)

This result holds also in the non-perturbative strong-driving regime, as long
as $T_{\phi}$ does not exceed $1/W$. Once that happens, $R$ is no longer well-
defined. In this regime $T_{\phi}$ can continue to increase, but that increase
is no longer due to the movement of weight to frequencies greater than $W$.

##  III Quantum Model for noise sources

Since the origin of $1/f$ noise is not well known in most cases, we need to
model the TLS as generally as possible. Here we introduce a quantum model of a
single TLS and compute its consequences for $T_{\phi}$.

It is particularly important to include the possibility that the TLS is a
quantum object in the sense that coherent superpositions of the two levels are
allowed, unlike the model of Sec. II. The key question is when the classical
SR behavior of the TLS is recovered when different parameters that govern the
quantumness of the model are varied.

###  III.1 Lindblad Equation

A single TLS is an open quantum system coupled to phonons (and possibly other
baths). The dynamics are described by the Lindblad equation [27]:

$\frac{d\rho}{dt}=-i[H,\rho]+\sum_{j}\left(L_{j}\rho
L_{j}^{\dagger}-\frac{1}{2}\left\\{L_{j}^{\dagger}L_{j},\rho\right\\}\right)$
(11)

where $\rho$ is the $2\times 2$ density matrix for the TLS and $\hbar=1$.

The off-diagonal elements $\rho_{01}$ and $\rho_{10}$ of the density matrix
have no classical analog.

The Hamiltonian is

$H=\frac{1}{2}(\epsilon\sigma_{z}+\Delta\sigma_{x})+(\alpha_{z}\sigma_{z}+\alpha_{x}\sigma_{x})\cos\omega_{d}t$
(12)

where $\epsilon$ is the energy difference between the two levels, $\Delta$ is
the tunneling amplitude, $\alpha_{z}$ and $\alpha_{x}$ are the driving
strengths that modulate the energy difference and the tunneling, respectively.
In practice, $\alpha_{z}$ and $\alpha_{x}$ are couplings that depend on the
properties of the qubit as well as the intensity of the driving field.
$\sigma_{i}$ are the Pauli operators and $\omega_{d}$ is the driving
frequency.

The Lindblad operators are

$L_{1}=\sqrt{\gamma}\begin{pmatrix}0&1\\\
0&0\end{pmatrix},\,L_{2}=\sqrt{\kappa}\begin{pmatrix}0&0\\\
1&0\end{pmatrix},\,L_{3}=\sqrt{\eta}\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}$
(13)

where $\gamma$, $\kappa$, and $\eta$ denote the relaxation, excitation, and
pure dephasing rates, respectively. If $\epsilon<0$ then $|0\rangle$ and
$|1\rangle$ are the ground and excited states respectively.

We now change to variables whose physical meaning is more transparent. As in
Sec. II, we use $s_{z}(t)=\rho_{00}(t)-\rho_{11}(t)$, and it is proportional
to the electric dipole moment in the case of charged fluctuators that produce
charge noise. The off-diagonal asymmetry is $p(t)=\rho_{01}(t)-\rho_{10}(t)$
and $q(t)=\rho_{01}(t)+\rho_{10}(t)$ is a phase measure. We also define the
total relaxation rate: $\Gamma=(\gamma+\kappa)/2$, and the relative rate:
$\lambda=\gamma-\kappa$.

The detailed balance condition is $2\Gamma\tanh(\epsilon/2k_{B}T)=\lambda$,
and we enforce this below in all our calculations, choosing $T=0.1$K.

Changing to these variables, we get the following coupled differential
equations from the Lindblad equation:

$\displaystyle\frac{ds_{z}(t)}{dt}$ $\displaystyle=-2\Gamma
s_{z}(t)+\lambda-(\Delta+2\alpha_{x}\cos\omega_{d}t)p(t)$ (14)
$\displaystyle\frac{dp(t)}{dt}$
$\displaystyle=(\Delta+2\alpha_{x}\cos\omega_{d}t)s_{z}(t)$
$\displaystyle\quad-(\Gamma+2\eta)p(t)-(\epsilon+2\alpha_{z}\cos\omega_{d}t)q(t)$
$\displaystyle\frac{dq(t)}{dt}$
$\displaystyle=(\epsilon+2\alpha_{z}\cos\omega_{d}t)p(t)-(\Gamma+2\eta)q(t).$

###  III.2 Classical Limit

Comparing Eqs. 14 to Eq. 1 in the absence of driving we see that with the
identifications $2\Gamma=W$ and $\lambda=\delta W$ that the Lindblad equations
we use reduce to the classical case when $p(t)=q(t)=0$ and $W_{+}$ and $W_{-}$
are independent of time. Furthermore, even when $p(t)$ and $q(t)$ are
included, if $\Gamma+2\eta$ is larger than the other parameters in the
problem, both $p(t)$ and $q(t)$ decay quickly and again the Lindblad equation
reduces to the classical limit after a short time (which is $T_{2}$).

While large $\Gamma+2\eta$ leads to the classical limit, large $\epsilon$ and
$\Delta$ increase the quantumness since they give rise to the Larmor
precession about the axis in the
$(x,y,z)=(\Delta,0,\epsilon)/\sqrt{\epsilon^{2}+\Delta^{2}}$. This aspect of
the quantum dynamics is absent in Eq. 1, a classical equation. The Larmor
frequency $\sqrt{\epsilon^{2}+\Delta^{2}}$ is an additional frequency
parameter in the problem, and thus we expect the PSD to have a complicated
frequency structure, in contrast to the classical case, where only the
harmonics of $\omega_{d}$ are involved. We will investigate this in Sec. V.

In addition, there are quantum effects in the structure of the Lindblad
equations themselves in that the driving terms are not diagonal in the
equation for $s_{z}$. For driving to exist, $p(t)$ must be present. This in
turn implies the existence of off-diagonal elements in the density matrix,
which means quantum coherence is present.

No closed-form solutions are known to exist for Eqs. 14, so we solve them
numerically. Ref. [28] treats some analytic methods that describe the dynamics
of a TLS under a strong driving field, but in a somewhat simpler model.

###  III.3 Computational Details

The focus is on $s_{z}(t)$ because it is proportional to the dipole moment
produced by the TLS. Since $s_{z}(t)=\langle\sigma_{z}\rangle=$
Tr$[\rho(t)\sigma_{z}]$ is an expectation value, the correlation function for
time duration $\tau$ can be regarded as $\pm s_{z}(\tau)$. It does depend on
the initial state, $|0\rangle$ or $|1\rangle$, of an ensemble. We are mainly
going to deal with $\epsilon$ much smaller than $k_{B}T$, and we find that in
the numerical calculations, the memory is sufficiently short that we can fix
the initial state to $|0\rangle$. (For further details of the computations,
see the Supplemental Material). Thus we can determine the quantity relevant
for decoherence, which is $C_{0}\int s_{z}(\tau)e^{-i\omega\tau}d\tau$.
Therefore, we will look into the changes in $S(\omega)$ and $T_{\phi}$ caused
by the driving field. To compute $S(\omega)$ we solve Eq. 14 numerically and
then use the fast Fourier transform on $s_{z}(t)$.

Since all of the interesting parameters may be expressed as frequencies when
$\hbar=k_{B}=1$, we set $\Gamma=1$ so that the basic frequency is the
switching rate of the TLS. This fixes the choice of units for all of the
quantities that appear in Eqs. 14 as well as all other frequencies mentioned
below.

$s_{z}(t)$ is sampled with frequency $f_{s}=10^{4}$ for a total simulation
time $T_{m}=10^{2}$ ($10^{6}$ points) for the single TLS calculations and
$f_{s}=10^{3}$ for $T_{m}=10^{3}$ (again $10^{6}$ points) for multiple TLSs.
This is sufficient for our purposes since we will only look at the changes in
$S(\omega)$ for $\omega<f_{s}/2$.

##  IV Effect of the TLS dephasing time

The basic ideas of this work are already evident at the single TLS level. We
assume in this section that the energy levels of the TLS are degenerate. This
allows a simple comparison with classical SR since the existence of that
phenomenon does not depend on the finiteness of the energy level separation.
In the quantum case, that separation leads to oscillatory behavior that will
be investigated in Sec. V.

###  IV.1 Parameters

We begin with the dephasing parameter $\eta=1$, which means that dephasing is
appreciable but not too strong. (Recall that the units are set by $\Gamma=1$)
This choice is suggested by studies on high-frequency TLSs having similar or a
few orders of magnitude smaller $\eta$ than $\Gamma$ [29, 30]. $k_{B}T/\hbar$
is normally much greater than the typical switching rates of TLS, and when we
apply detailed balance we find $\lambda=0$. In this section, we also take
$\Delta=0$, and $\epsilon=0$. The two energy levels of the TLS are degenerate
so there is no Larmor precession. Precession effects will be treated in the
next section. Here the quantum effects arise only from the mixing of
$s_{z}(t)$ with $p(t)$ and $q(t)$ as explained above. We will vary the driving
frequency $\omega_{d}$ and overall driving strength but enforce the
relationship $\alpha_{x}=0.5\alpha_{z}$. We may think of this relationship as
off-axis driving since the precession axis of the driving field does not align
with the $z$-axis that is associated with measurable quantities.

###  IV.2 Power spectral density

In general, we find that $\alpha_{x}$ is much more important than
$\alpha_{z}$; there will be no driving effect if $\alpha_{x}\ll\alpha_{z}$.
For further discussion of their relative magnitudes, see the Supplemental
Material.

Figure 2: (a) Shift of the PSD for a single TLS varying $\alpha_{z}$ with
fixed $\omega_{d}=10^{1}$. These PSDs correspond to the three boxed simulation
points along the vertical axis in Fig. 3. (b) Shift of PSD for a single TLS
varying $\omega_{d}$ with fixed $\alpha_{z}=10^{4}$. These PSDs are offset by
2 orders of magnitude for clarity and correspond to the three boxed simulation
points along the horizontal axis in Fig. 3. Color-coded dashed lines mean the
offset power at $\omega=10^{-2}$ when the driving field is not applied. All
frequencies are in units of $\Gamma$, the width of the Lorentzian centered at
$\omega=0$. Figure 3: Increase in dephasing time of a qubit when a single TLS
is driven by the oscillating electric field with strength $\alpha_{z}$ and
frequency $\omega_{d}$. The dephasing time with driving field ($T_{\phi,d}$)
is divided by that without driving field ($T_{\phi,0}$). Only 49 values of
$T_{\phi,d}/T_{\phi,0}$ with the combination of 7 different strengths
($10^{-2}$ to $10^{4}$) and frequencies ($10^{-2}$ to $10^{4}$) are shown.
White dashed boxes are for parameters used for the curves in Fig. 2.

Fig. 2 gives the main results for the PSD.

Fig. 2(a) shows the PSD when varying the driving strength $\alpha_{z}$. The
first important point is that the strength of the central Lorentzian does
indeed decrease as $\alpha_{z}$ increases. This means that this basic feature
of classical SR is preserved in the quantum density matrix formalism. The
decrease is appreciable, quite visible on the log-log plot in the frequency
range $\omega<1$. However, there are some differences between the results and
what might be expected from classical SR. We do not see a sharp peak at
$\omega_{d}=10$. Instead, the weight for moderate driving is spread across a
range of frequencies near $\omega_{d}$. For strong driving, the weight moves
up still further in frequency.

Fig. 2(b) shows the PSD when the driving strength is strong
($\alpha_{z}=10^{4}$) and fixed while $\omega_{d}$ is varied. The curves are
offset, so the reduction in weight of the Lorentzian at low frequency is not
easy to see. The weight at high frequencies moves to the right as $\omega_{d}$
is increased. This is in qualitative agreement with classical SR but again
this weight is spread across a range of frequencies: there is no evidence of
isolated sharp peaks. The structure at the far right in the small $\omega_{d}$
case is due to sampling error. Overall, the redistribution ratio $R$ is a
function of both the driving strength and the driving frequency.

In practice, the decrease of $R$ may be mainly limited by practical
considerations. The frequency dependence is more interesting and shows the
same trade-off mentioned in Sec. II. At first, the shift of the weight ($R$
decreases) increases with frequency as the TLS attempts to follow the driving
frequency, in somewhat the same fashion as a driven harmonic oscillator.
However, as the driving frequency continues to increase, the driving force
averages out and the effect decreases. We conclude that in the case where the
Larmor precession is absent ($\epsilon=\Delta=0$) and intrinsic dephasing is
moderate ($\eta=1$), quantum effects do not destroy the phenomenon of SR. The
main effect is to greatly complicate the frequency structure of the PSD at
high frequencies. The single delta-function peak at $\omega_{d}$ converts to
comb-like structures at higher frequencies.

###  IV.3 Qubit dephasing time

Any change in $S(\omega)$ will be reflected in the dephasing time $T_{\phi}$
of nearby qubits. We use Eq. 8 with a Gaussian window function to compute
$T_{\phi}$ as a function of driving strength and frequency. To see the effects
of driving, we compute the ratio of $T_{\phi,d}$, the dephasing time with
driving, to $T_{\phi,0}$, that in the absence of driving. The overall coupling
strength $C_{0}$ is not so important here since we focus on the relative
suppression of decoherence. However, to give an idea of experimentally
relevant parameters in dimensionful units, if we take a coupling strength
$C_{0}=10^{12}$ Hz2 and a width of $W_{0}=1$ Hz, we will find
$T_{\phi,0}=1.00$ $\upmu$s, a typical value in spin qubit experiments. We
emphasize once more that these are rough estimates since Eq. 8 is an
approximation.

In the simulations in Fig. 3, $C_{0}=1$ is used, resulting in
$T_{\phi,0}=2.31$. The key point is that to have a significant increase in
$T_{\phi}$, it is critical that $\omega_{d}>1/T_{\phi,0}$. If $\omega_{d}$ is
smaller than this, the weight in the PSD may be shifted, but the shifted
weight still contributes to the decoherence of the qubit, which is not what we
want. Very large values of $\omega_{d}$ are also not favorable. All of this is
in agreement with classical SR.

As shown in Fig. 3, the increase in $T_{\phi}$ starts when
$\alpha_{z}>\omega_{d}>\Gamma$, enhancing $T_{\phi}$ by a factor of up to
1.79. There is no improvement in $T_{\phi}$ when $\alpha_{z}<\Gamma$ or
$\alpha_{z}<\omega_{d}$. This is in accord with the changes in $S(\omega)$
discussed above: the time scale of the driving strength should be shorter than
the average relaxation time of the TLS to see the maximum effect.

###  IV.4 Summary

Fig. 4 shows the results of varying the parameter $\eta$, which controls the
dephasing, that is the decay of the off-diagonal parts of the density matrix
$p$ and $q$. Note that every curve crosses every other one: this is a
consequence of the preservation of the total weight.

$\eta$ can be thought of as a purely quantum parameter that complements
$\Gamma$, which has a classical analog, as noted above. In the $\hat{z}$
basis, both contribute to $T_{2}$ processes, while $\Gamma$ alone contributes
to $T_{1}$ processes.

When $\eta$ is relatively small ($\eta/\omega_{d}\leq 1$) we see the effects
of classical SR, as already noted. Roughly speaking, the shift in weight can
be followed by looking at the asymptotic value of the PSD at small $f$. As
$\eta$ increases, this value first decreases as the width of the Low-frequency
Lorentzian increases, and the height decreases. The concomitant noisy behavior
at higher frequencies is pronounced. As $\eta$ increases still further, past
$10\,\omega_{d}$, decay of $p$ and $q$ is so rapid that no oscillations from
the driving are possible. Exponential decay at a rate $\propto\Gamma$ results
and there is no shift in weight. The net effect is a non-monotonic shift in
weight as a function of $\eta$, with the maximum around
$\eta/\omega_{d}\approx 10$. The black arrows in Fig. 4 indicate this. This
non-monotonicity is a purely quantum effect. We expect that it would manifest
itself in an experiment as a non-monotonic behavior in $T_{\phi,d}/T_{\phi,0}$
as a function of the controllable parameter $\omega_{d}$.

Figure 4: Power spectral density for different values of $\eta$, the dephasing
parameter of the TLS. This is the case that the energy levels of the TLS are
degenerate: $\epsilon=0$. The driving strengths are
$\alpha_{z}=\alpha_{x}=100$, and the driving frequency is $\omega_{d}=10$. The
black arrows go in the direction of increasing $\eta$ and demonstrate the non-
monotonic dependence of the redistribution ratio $R$ on the dephasing rate of
the TLS.

Figure 5: Enhancement of the dephasing time as a function of driving strength
$\alpha_{z}$ and driving frequency $\omega_{d}$ of a qubit with energy level
differences (a) $\epsilon=0.1$ and (b) $\epsilon=10$. $\Gamma=1$, $\eta=1$,
$\alpha_{x}=0.5\alpha_{z}$. The slow precession in (a) does not change the
overall pattern of the enhancement from the classical case. The fast-
precessing system in (b) shows a new peak in the enhancement when the
precession frequency $\epsilon$ matches the scaled driving amplitude
$\alpha_{z}/\Gamma$.

##  V Effect of Larmor precession of the TLS

We now introduce a nonzero energy level separation $\epsilon$ into the
Hamiltonian. The two levels of a classical TLS also have different energies in
general. However, in the quantum case, this difference also implies a Larmor
precession, an entirely nonclassical effect.

We only consider the special case that $\Delta=0$. This means that the
precession axis is the $z$-axis of the Bloch sphere. Since we are looking at
the correlation function of $s_{z}$, this means that the effect of the
precession is somewhat indirect, as a rotation about the $z$-axis does not
affect $s_{z}$. This implies that in the absence of driving there is no peak
in the PSD at $\Delta$.

In Fig. 5 we plot the $T_{\phi}$ enhancement as a function of the driving
strength and the driving frequency for $\epsilon=0.1$ (slow precession) and
$\epsilon=10$ (fast precession). Slow precession is not enough to change the
basic picture of essentially classical SR that we saw in the previous section.
The enhancement still depends monotonically on the driving strength and there
is the same peak in the enhancement as a function of driving frequency. When
the precession is fast the picture changes in one very important respect. A
peak in the enhancement as a function of driving strength appears at
$\alpha_{z}\approx\epsilon$.

This peak can serve as an experimental signature of the energy level spacing
made possible by a quantum effect.

##  VI PSD and dephasing time of multiple TLSs

Figure 6: Shift of the PSD of multiple TLSs as $\omega_{d}$ changes with fixed
$\alpha_{z}=10^{3}$. These PSDs correspond to the two simulation points in
Fig. 7, respectively. All frequencies are in units of $\Gamma_{mid}$, the
switching rate in the logarithmic center of the seven TLS that produce the
noise. (a) Here $\omega_{d}$ is in the range of the switching rates of the
TLS. (b) Here $\omega_{d}$ is above the range of the switching rates of the
TLS. Note the suppression of the low-frequency noise in panel (a). This is
nearly absent in panel (b). Figure 7: Increase in dephasing time of a qubit
when multiple TLSs are driven by the oscillating electric field with strength
$\alpha_{z}$ and frequency $\omega_{d}$. The dephasing time with driving field
($T_{\phi,d}$) is divided by that without driving field ($T_{\phi,0}$). 49
values of $T_{\phi,d}/T_{\phi,0}$ with the combination of 7 different
strengths and frequencies are shown. $\Gamma_{mid}$ is the frequency in the
logarithmic middle of the undriven PSD. White dashed boxes are the simulation
points used in Fig. 6.

The phenomenon of decoherence suppression has been established for the case of
a single TLS. In real systems many TLSs with a distribution of relaxation
times are acting on the qubit, giving rise to PSDs that are generally not
Lorentzian. This suggests the question of whether a single driving source with
a single driving frequency can still suppress qubit decoherence.

We do not attempt a complete analysis of this question in this paper, which
would involve optimizing the driving parameters for each type of PSD. Instead,
we merely try to settle the issue of whether NSD can mitigate the effects of
the best-known PSD: $1/f$ noise. This is produced by a set of TLS with a log-
uniform distribution of Lorentzian widths. Somewhat surprisingly, only one TLS
per frequency decade is required to give a 1/f spectrum quite accurately.
Motivated by this and computational convenience, we use 7 TLSs in this
section. We set $\Gamma=\eta$ for each TLS. Their switching rates are set as
$10^{-6}$, $10^{-5}$,$10^{-4}$,$10^{-3}$,$10^{-2}$,$10^{-1}$, and $10^{0}$.
Thus for this calculation, the unit of frequency is set equal to the fastest
switching rate. We also define $\Gamma_{mid}=10^{-3}$, which corresponds to
logarithmic middle of the switching rates. Other parameters are set to those
used in Sec. IV. Each TLS is simulated separately using the same method as
above, and the individual power spectra $S_{i}(\omega)$ are added to give the
total PSD, $S(\omega)=\sum_{i}S_{i}(\omega)$, which then yields the
fluctuating field at the qubit.

The results are given in Fig. 6, where we plot the total PSDs with fixed
driving strength and two different driving frequencies. We include comparison
plots with no driving. The latter gives $1/f$ noise; the slight curvature in
the plots is due to the fact that we only include 7 TLS, so at high
frequencies the PSD has a $1/f^{2}$ behavior. The changes in the PSD are more
subtle than in the single TLS case. The spectral weight can be moved to higher
frequencies, but only in a limited range of $\omega_{d}$. This is a
consequence of the fact that for the driving to have an effect on a TLS there
needs to be an approximate match between the switching rate of the TLS and the
driving frequency $\omega_{d}$. This cannot be achieved for all of the TLSs
when the noise is broadband since we limit ourselves to a single driving
frequency.

It is clearly more effective to set $\omega_{d}$ in the frequency range of the
TLSs, as in Fig.6 (a), where $\omega_{d}=10^{-1}$, rather than above that
range, as in Fig.6 (b), where $\omega_{d}=10^{2}$. The difference is most
clearly understood by looking at small $f/\Gamma_{mid}$, where the suppression
of the low-frequency noise noise is markedly suppressed in Fig. 6 (a) where
$\omega_{d}/\Gamma_{mid}=10^{2}$, but much less so in Fig. 6 (b), where
$\omega_{d}/\Gamma_{mid}=10^{5}$. This is entirely in agreement with the
above-mentioned importance of the matching of $\omega_{d}$ with switching
rates.

The results for the change in the dephasing time for 7 TLSs are given in Fig.
7. We see a qualitatively similar effect to that in Fig. 3, but find the
maximum value of $T_{\phi,d}/T_{\phi,0}$ to be about 2.43 when
$\alpha_{z}=10^{3}$ and $\omega_{d}=10^{-1}$. This is quite intriguing because
$\omega_{d}$ is smaller than $1/T_{\phi,0}=2.27$ and inside the noise
bandwidth ($10^{-6}-1$). This strongly suggests that an optimal driving
frequency should exist for a given $1/f$-like noise spectrum. It is also most
likely possible to do better using non-monochromatic driving. However, the
detailed optimization of the driving field is beyond the scope of this paper.

##  VII Discussion

Driving the TLS noise sources does indeed change the PSD by moving weight to
high frequencies where its effect on $T_{\phi}$ is reduced and qubit coherence
is thereby enhanced. This may be a promising way to improve coherence in a
system where dephasing from TLS is an important issue. The weight shift occurs
already at the classical level - it is present also in the theory of classical
SR. It does survive in the quantum case. This can happen when the TLS is
effectively classical because of its own rapid dephasing. This is in fact the
usual picture of a TLS, and in this paper we showed how this picture is
recovered in one limit of the quantum equations.

An important conclusion of the analysis was that the dependence of $T_{\phi}$
on the driving frequency is non-monotonic, with a peak in the enhancement when
$\omega_{d}$ is comparable to or somewhat less than the TLS dephasing rate
$\eta$. This provides a way for experiments to determine the rate of TLS
dephasing, a quantity that is otherwise difficult to determine.

Precession of the TLS Bloch vector is another quantum effect that can affect
the noise produced by the TLS. It comes from $\epsilon$, the difference in
energy of the two levels. When the precession is slow compared to the
relaxation rate $\Gamma$ or to the driving strength parameter $\alpha_{z}$
then it does not affect the enhancement of qubit coherence. However, when
$\epsilon$ exceeds $\Gamma$ then there is a resonance when
$\epsilon\approx\alpha_{z}$. Observing this resonance as $\alpha_{z}$ is
adjusted would give a measure of $\epsilon$. To see these effects on the PSD,
it would be necessary that only one TLS dominates the noise at the qubit,
which is of course not the typical situation. However, specially designed
experiments have succeeded in doing this [31].

One objection to NSD is that it may be difficult in practice to drive the TLS
electromagnetically without also driving the qubit at some level, with
possible deleterious effects. For example, NSD may heat the sample. This can
induce a shift of resonant frequency for a spin qubit [32, 33]. Sonic driving
is also a possibility, though also here one may anticipate unwanted qubit
driving, since elastic vibrations would, for example, modulate tunneling
matrix elements in two-qubit gates that depend on the exchange interaction.
However, the driving frequencies $\omega_{d}$ that are of interest are those
that tend to be in the range of the frequencies of the slow noise that produce
qubit dephasing, not the much higher frequencies characteristic of qubit
operating frequencies. This extreme mismatch implies that it should be
possible to minimize the direct effect of the driving field on the qubits. We
note also that dressed qubit protocols may be a useful tool for treating those
effects that still remain. [34, 35].

The physics of the TLS, the qubits, and their couplings clearly involve many
parameters, and this paper covers only a small portion of the whole space. For
example, there is a regime where the qubit-TLS coupling $C_{0}$, in our
notation, is large enough that the qubit can affect the dynamics of the TLS.
Such effects have been observed experimentally [36]. The relatively simple
equation we used for $T_{\phi}$ may not be valid in this strong coupling
regime [37]. We considered only one off-axis direction for the driving field,
and there may be interesting physics that was missed because of this
limitation. Temperature was only included by enforcing detailed balance
between the upward and downward transitions. This problem can certainly
benefit from further investigation.

##  VIII Acknowledgments

We thank Leah Tom for providing COMSOL simulation results in the Supplemental
Material. This research was partly conducted while visiting the Okinawa
Institute of Science and Technology (OIST) through the Theoretical Sciences
Visiting Program (TSVP). This research was sponsored by the Army Research
Office (ARO) under Awards No. W911NF-17-1-0274 and No.
W911NF-W911NF-23-1-0110. The views, conclusions, and recommendations contained
in this document are those of the authors and are not necessarily endorsed nor
should they be interpreted as representing the official policies, either
expressed or implied, of the Army Research Office (ARO) or the U.S.
Government. The U.S. Government is authorized to reproduce and distribute
reprints for Government purposes notwithstanding any copyright notation
herein. YC acknowledges support from the Army Research Office through Grant
No. W911NF-23-1-0115.

## References

  * Paladino  _et al._ [2014] E. Paladino, Y. Galperin, G. Falci, and B. Altshuler, 1/f noise: Implications for solid-state quantum information, Reviews of Modern Physics 86, 361 (2014).
  * Müller  _et al._ [2019] C. Müller, J. H. Cole, and J. Lisenfeld, Towards understanding two-level-systems in amorphous solids: insights from quantum circuits, Reports on Progress in Physics 82, 124501 (2019).
  * Connors  _et al._ [2019] E. J. Connors, J. Nelson, H. Qiao, L. F. Edge, and J. M. Nichol, Low-frequency charge noise in Si/SiGe quantum dots, Phys. Rev. B 100, 165305 (2019).
  * Freeman  _et al._ [2016] B. M. Freeman, J. S. Schoenfield, and H. Jiang, Comparison of low frequency charge noise in identically patterned Si/SiO2 and Si/SiGe quantum dots, Applied Physics Letters 108, 253108 (2016).
  * Zimmerman  _et al._ [2014] N. M. Zimmerman, C.-H. Yang, N. S. Lai, W. H. Lim, and A. S. Dzurak, Charge offset stability in Si single electron devices with Al gates, Nanotech. 25, 405201 (2014).
  * Rudolph  _et al._ [2019] M. Rudolph, B. Sarabi, R. Murray, M. Carroll, and N. M. Zimmerman, Long-term drift of Si-MOS quantum dots with intentional donor implants, Scientific reports 9, 7656 (2019).
  * Takeda  _et al._ [2013] K. Takeda, T. Obata, Y. Fukuoka, W. M. Akhtar, J. Kamioka, T. Kodera, S. Oda, and S. Tarucha, Characterization and suppression of low-frequency noise in Si/SiGe quantum point contacts and quantum dots, Applied Physics Letters 102, 123113 (2013).
  * Suter and Álvarez [2016] D. Suter and G. A. Álvarez, Colloquium: Protecting quantum information against environmental noise, Reviews of Modern Physics 88, 041001 (2016).
  * Szańkowski  _et al._ [2017] P. Szańkowski, G. Ramon, J. Krzywda, D. Kwiatkowski, and Ł. Cywiński, Environmental noise spectroscopy with qubits subjected to dynamical decoupling, Journal Of Physics-Condensed Matter 29, 333001 (2017).
  * De Lange  _et al._ [2012] G. De Lange, T. Van Der Sar, M. Blok, Z.-H. Wang, V. Dobrovitski, and R. Hanson, Controlling the quantum dynamics of a mesoscopic spin bath in diamond, Scientific Reports 2, 1 (2012).
  * Knowles  _et al._ [2014] H. S. Knowles, D. M. Kara, and M. Atatüre, Observing bulk diamond spin coherence in high-purity nanodiamonds, Nature Materials 13, 21 (2014).
  * Bauch  _et al._ [2018] E. Bauch, C. A. Hart, J. M. Schloss, M. J. Turner, J. F. Barry, P. Kehayias, S. Singh, and R. L. Walsworth, Ultralong dephasing times in solid-state spin ensembles via quantum control, Physical Review X 8, 031025 (2018).
  * Bluvstein  _et al._ [2019] D. Bluvstein, Z. Zhang, C. A. McLellan, N. R. Williams, and A. C. B. Jayich, Extending the quantum coherence of a near-surface qubit by coherently driving the paramagnetic surface environment, Physical Review Letters 123, 146804 (2019).
  * Joos  _et al._ [2022] M. Joos, D. Bluvstein, Y. Lyu, D. Weld, and A. Bleszynski Jayich, Protecting qubit coherence by spectrally engineered driving of the spin environment, npj Quantum Information 8, 47 (2022).
  * You  _et al._ [2022] X. You, Z. Huang, U. Alyanak, A. Romanenko, A. Grassellino, and S. Zhu, Stabilizing and improving qubit coherence by engineering the noise spectrum of two-level systems, Physical Review Applied 18, 044026 (2022).
  * Lisenfeld  _et al._ [2023] J. Lisenfeld, A. Bilmes, and A. V. Ustinov, Enhancing the coherence of superconducting quantum bits with electric fields, npj Quantum Information 9, 8 (2023).
  * Gammaitoni  _et al._ [1998] L. Gammaitoni, P. Hänggi, P. Jung, and F. Marchesoni, Stochastic resonance, Reviews of Modern Physics 70, 223 (1998).
  * J.D.  _et al._ [2017] B. J.D., B. A., W. G., U. A.V., and L. J, Transmission-line resonators for the study of individual two-level tunneling systems, Appl. Phys. Lett. 111, 112601 (2017).
  * McNamara and Wiesenfeld [1989] B. McNamara and K. Wiesenfeld, Theory of stochastic resonance, Physical Review A 39, 4854 (1989).
  * Shneidman  _et al._ [1994] V. A. Shneidman, P. Jung, and P. Hanggi, Weak-noise limit of stochastic resonance, Physical Review Letters 72, 2682 (1994).
  * Gammaitoni  _et al._ [1995] L. Gammaitoni, F. Marchesoni, E. Menichella-Saetta, and S. Santucci, Stochastic resonance in the strong-forcing limit, Physical Review E 51, R3799 (1995).
  * Löfstedt and Coppersmith [1994a] R. Löfstedt and S. Coppersmith, Stochastic resonance: Nonperturbative calculation of power spectra and residence-time distributions, Physical Review E 49, 4821 (1994a).
  * Wagner  _et al._ [2019] T. Wagner, P. Talkner, J. C. Bayer, E. P. Rugeramigabo, P. Hänggi, and R. J. Haug, Quantum stochastic resonance in an ac-driven single-electron quantum dot, Nature Physics 15, 330 (2019).
  * Hänze  _et al._ [2021] M. Hänze, G. McMurtrie, S. Baumann, L. Malavolti, S. N. Coppersmith, and S. Loth, Quantum stochastic resonance of individual Fe atoms, Science Advances 7, eabg2616 (2021).
  * Löfstedt and Coppersmith [1994b] R. Löfstedt and S. Coppersmith, Quantum stochastic resonance, Physical Review Letters 72, 1947 (1994b).
  * Grifoni and Hänggi [1996] M. Grifoni and P. Hänggi, Nonlinear quantum stochastic resonance, Physical Review E 54, 1390 (1996).
  * Manzano [2020] D. Manzano, A short introduction to the Lindblad master equation, AIP advances 10 (2020).
  * Saiko  _et al._ [2016] A. Saiko, S. Markevich, and R. Fedaruk, Dissipative two-level systems under ultrastrong off-resonant driving, Physical Review A 93, 063834 (2016).
  * Lisenfeld  _et al._ [2010] J. Lisenfeld, C. Müller, J. H. Cole, P. Bushev, A. Lukashenko, A. Shnirman, and A. V. Ustinov, Measuring the temperature dependence of individual two-level systems by direct coherent control, Physical Review Letters 105, 230504 (2010).
  * Lisenfeld  _et al._ [2016] J. Lisenfeld, A. Bilmes, S. Matityahu, S. Zanker, M. Marthaler, M. Schechter, G. Schön, A. Shnirman, G. Weiss, and A. V. Ustinov, Decoherence spectroscopy with individual two-level tunneling defects, Scientific Reports 6, 1 (2016).
  * Ye  _et al._ [2024] F. Ye, A. Ellaboudy, D. Albrecht, R. Vudatha, N. T. Jacobson, and J. M. Nichol, Characterization of individual charge fluctuators in Si/SiGe quantum dots, arXiv preprint arXiv:2401.14541 (2024).
  * Undseth  _et al._ [2023] B. Undseth, O. Pietx-Casas, E. Raymenants, M. Mehmandoost, M. T. Mądzik, S. G. Philips, S. L. De Snoo, D. J. Michalak, S. V. Amitonov, L. Tryputen, _et al._ , Hotter is easier: unexpected temperature dependence of spin qubit frequencies, Physical Review X 13, 041015 (2023).
  * Choi and Joynt [2024] Y. Choi and R. Joynt, Interacting random-field dipole defect model for heating in semiconductor-based qubit devices, Physical Review Research 6, 013168 (2024).
  * Seedhouse  _et al._ [2021] A. E. Seedhouse, I. Hansen, A. Laucht, C. H. Yang, A. S. Dzurak, and A. Saraiva, Quantum computation protocol for dressed spins in a global field, Physical Review B 104, 235411 (2021).
  * Hansen  _et al._ [2021] I. Hansen, A. E. Seedhouse, A. Saraiva, A. Laucht, A. S. Dzurak, and C. H. Yang, Pulse engineering of a global field for robust and universal quantum computation, Physical Review A 104, 062415 (2021).
  * Brox  _et al._ [2011] H. Brox, J. Bergli, and Y. M. Galperin, Effects of external driving on the coherence time of a josephson junction qubit in a bath of two-level fluctuators, Physical Review B 84, 245416 (2011).
  * Bergli  _et al._ [2009] J. Bergli, Y. M. Galperin, and B. Altshuler, Decoherence in qubits due to low-frequency noise, New Journal of Physics 11, 025002 (2009).


# Nonlinear clogging of a rectangular slit by a spherical soft particle

Charles Paul Moore  Institut Pasteur, Université Paris Cité, Physical
microfluidics and Bioengineering, 25-28 Rue du Dr Roux, 75015 Paris, France
LadHyX, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, 91120
Palaiseau, France  Julien Husson  LadHyX, CNRS, Ecole Polytechnique, Institut
Polytechnique de Paris, 91120 Palaiseau, France  Arezki Boudaoud  LadHyX,
CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau,
France  Gabriel Amselem  LadHyX, CNRS, Ecole Polytechnique, Institut
Polytechnique de Paris, 91120 Palaiseau, France  Charles N. Baroud  Institut
Pasteur, Université Paris Cité, Physical microfluidics and Bioengineering,
25-28 Rue du Dr Roux, 75015 Paris, France  LadHyX, CNRS, Ecole Polytechnique,
Institut Polytechnique de Paris, 91120 Palaiseau, France

###### Abstract

The capture of a soft spherical particle by a rectangular slit leads to a non-
monotonic pressure-flow rate relation at low Reynolds number. In the presence
of the trapped particle the flow-induced deformations focus the streamlines
and pressure drop to a small region. This increases the resistance to flow by
several orders of magnitude as the driving pressure is increased. As a result
two regimes are observed: a flow-dominated regime for small particle
deformations, where flow rate increases with pressure, and an elastic-
dominated regime in which solid deformations block the flow.

Solid deformations can have a major influence on fluid flow in the regime of
strong fluid-structure coupling. A wide range of complex behaviors, including
oscillations that lead to catastrophic failures or blockages, can be observed
when the inertial effects of the fluid and solid couple together [1]. Recently
an increased attention has been paid to the effet of a low Reynolds number
flow on an elastic fiber [2], with particular focus on the transport [3] and
deformation [4] of slender elastic fibers by the flow [5].

Two-way flow-structure coupling has been studied in the case of deformable
tubes, in which elastic deformations of the tube walls led to strong
modifications of the fluid flow. Wall elasticity effects were shown to
stabilize the flow distribution into a bifurcation [6] and to suppress the
emergence of viscous fingering in a Hele-Shaw cell [7]. More extreme cases
emerged in the case of air-liquid flows within flexible tubes, where surface
tension led to a complete airway closure [8]. These effects of wall elasticity
of the fluid flow have been used to explain sap flow in green plants [9] or to
create soft valves for technological applications [10, 11, 12].

A different class of problems for which the two-way coupling can lead to
extreme modifications of the flow consists of the case of a soft particle
being pushed into an orifice. This problem is encountered in many microfluidic
applications, such as during the flow and encapsulation of hydrogels [13, 14]
or for the characterization of cells and other soft materials [15, 16, 17, 18,
19, 20]. It is also closely related to clogging or sieving particles [21, 22,
23, 24, 25, 26, 27]. In many of these cases a spherical particle is pushed
into a rectangular slit, which allows for a leakage flow around the particle
[15, 17, 28, 19]. But in spite of the interest for applications, the physics
that determines the equilibrium between leaky flow and particle deformation
has not been explored. Here we describe this nonlinear relationship and
explain the low and high deformation regimes that are encountered.

Figure 1: (a) Sketch of the microfluidic device. The flow rate is imposed at
the inlet and divides into a thrupass and bypass channels. The micrograph
shows a gel bead trapped in a slit of width $w=25$ µm. A resistance diagram of
the microfluidic channel is included. (b) The flow rate $Q_{\rm thru}$ past
the microgel increases to a maximum and then decreases to zero as the inlet
flow rate $Q$ is increased. Each line represents a different trapped microgel.
Young’s modulus $E=2.1\pm 0.9\;\rm kPa$. Gel diameter: $d=83\pm 3\;\mu\rm m$.
Slit width: $w=25\pm 5\;\mu\rm m$. Channel height: $h=85\;\mu\rm m$. Inset:
reconstructed streamlines showing the flow division between the thrupass and
bypass channels. The value of $y/y_{\rm max}$ provides a measure of the flow
rates $Q_{\rm thru}$ and $Q_{\rm byp}$.

The experimental setup consisted in a microfluidic device with two parallel
channels: a thrupass line with a narrow slit to trap microgel beads, and a
bypass line of width $y_{\rm max}$, which provided a known hydraulic
resistance, see Fig. 1. The flow rate $Q$ at the inlet was imposed with a
syringe pump (Nemesys, Cetoni GMBH, Germany), and streamlines of the flow in
the device were reconstructed by tracking fluorescent particles (0.5 $\mu\rm
m$ Polybead, Polysciences Inc., Washington Pa) with PIVLab [29, 30], see inset
in Fig. 1b. The flow rate $Q_{\rm thru}$ in the thrupass channel was obtained
by measuring the position $y$ of the separatrix streamline, shown in the inset
of Fig. 1b. The separatrix divides the fluid going to the thrupass and the
bypass channels, according to: $Q_{\rm
thru}={\int_{0}^{y}U(y)dy}/{\int_{0}^{y_{\rm max}}U(y)dy}Q$, where $U(y)$ is
the theoretical flow speed in a rectangular channel of a height $h$ and width
$y_{\rm max}$ [31].

In the absence of a gel bead, increasing the flow rate $Q$ at the inlet led to
an increase of the flow rate $Q_{\rm thru}$, with the ratio between $Q_{\rm
thru}$ and $Q$ being given by the ratio of the hydrodynamic resistances
between the thrupass and bypass channels. When a gel bead was trapped by the
narrow slit, increasing the flow rate $Q$ at the inlet led to a non-monotonic
behavior of the flow rate $Q_{\rm thru}$ through the thrupass channel. At
first, increasing $Q$ led to an increase in $Q_{\rm thru}$, until a maximum
value of $Q_{\rm thru}$ was reached; increasing $Q$ beyond a critical value
then led to a decrease in the flow rate $Q_{\rm thru}$, see Fig. 1b. At high
enough values of $Q$, the gel bead plugged the thrupass channel entirely and
all the flow was redirected to the bypass channel. This non-monotonic
relationship between $Q_{\rm thru}$ and $Q$ reveals a strong non-linear
hydrodynamic resistance added by the gel bead.

Figure 2:  (a) Hydrodynamic resistance due to the gel $R_{\rm gel}$, as a
function of the imposed flow rate $Q$. Note that the resistance increases
faster for softer gels than stiffer gels. (b) Microgel displacement as a
function of the applied pressure. Each curve represents data collected on a
single microgel.

Experiments were repeated for poly(ethylene glycol) (PEG) beads of diameters
$d=80\;-145\;\mu\rm m$. The beads were made using droplet flow lithography by
adding PEG-diacrylate (PEG-DA) to a PEG solution and photo-polymerizing it in
a first microfluidic device [32]. The beads were stored off-chip until they
were re-injected into the current device. The microgel stiffness was
controlled by modulating the ratio of PEG to PEG-DA. The values of the
equivalent Young’s moduli were obtained using a microindentation technique
[33] and spanned two orders of magnitude: $E^{\star}=100\;\rm Pa$ to $21\;\rm
kPa$, where $E^{\star}=\frac{E}{{1-\nu^{2}}}$, $E$ is the conventional Young’s
modulus and $\nu$ is the Poisson ratio. This modulus is measured directly by
microindentation and accounts for contact loading rather than simple
compression/extension. Three different microfluidic traps were used, ranging
within $w=15-45\;\mu\rm m$, and microchannel heights $h$ chosen to
approximately match microgel diameters: $h=d\pm 10\%$.

At low Reynolds number, the pressure $\Delta P$ across a channel is
proportional to the flow $Q$ in the channel and the hydrodynamic resistance
$R$ of the channel: $\Delta P=QR$ [34]. Calling $R_{\rm thru}$ and $R_{\rm
byp}$ the known hydrodynamic resistances of the thrupass and bypass channels
in the absence of gel and $R_{\rm gel}$ the added hydrodynamic resistance of
the gel, we therefore have: $\frac{Q_{\rm byp}}{Q_{\rm thru}}=\frac{(R_{\rm
thru}+R_{\rm gel})}{R_{\rm byp}}$, which provides a way to compute the
resistance $R_{\rm gel}$ added by the gel, as well as the pressure drop across
the gel $P_{\rm gel}=R_{\rm gel}Q_{\rm thru}$ (see Fig. 1a).

The value of $R_{\rm gel}$ increased dramatically with flow rate, spanning
nearly 6 orders of magnitude when the flow rates covered 2 orders of magnitude
(see Fig. 2a). The rate of this increase depended on the gel elasticity:
stiffer gels ($E^{\star}\approx 10^{4}\;\rm Pa$) led to a slower increase in
resistance than softer gels ($E^{\star}\approx 10^{2}\;\rm Pa$). The increase
in resistance was associated with a displacement of the microgel as it
deformed and penetrated into the slit. This displacement was quantified by
measuring the distance $l$ traveled by the back of the gel, i.e. upstream of
the slit, with respect to its position in the absence of flow (see inset in
Fig. 2b). The value of $l$ increased with the pressure drop across the gel
$P_{\rm gel}$, see Fig. 2b. For a given pressure $P_{\rm gel}$ across the gel
bead, the softer gels penetrated more into the slit, and even more so when the
gap width was larger, see Fig. 2b.

To understand how the flow forces the hydrogel to deform, plug the slit and
increase the hydrodynamic resistance, simulations of a soft particle deforming
into a slit were performed on a quarter-setup using the software Abaqus. The
solid deformations were modeled by simulating an initially spherical gel bead
that was subjected to a negative pressure inside the slit. The resulting
deformed gel geometry (Fig. 3a) was exported to the software COMSOL
Multiphysics and the flow field in the channel containing the deformed gel was
simulated, providing the pressure distribution in the fluid everywhere in the
device. The uncoupled solid and fluid simulations were then iterated once by
updating the gel shape and then the pressure field.

Figure 3:  (a) A quarter section of the deformed microgel and trap for an
applied pressure $P=4\;\rm kPa$. Color represents relative nodal displacement
in the microgel. (b) Simulation of the flow around an undeformed gel (left),
and a gel deformed under a simulated pressure of $P_{\rm gel}=4\;\rm kPa$
(right). Color indicates the pressure normalized by the inlet pressure $P$.
Isobars are shown in white. Note how the deformed microgel focuses the
pressure drop into a narrow area. $x$ and $y$ directions represent the imaging
plane. (c) The measured gel hydrodynamic resistance $R_{\rm gel}$ as a
function of the area $\Lambda$ of the triangle through which the flow must
pass (see inset). $R_{\rm gel}$ and $\Lambda$ are respectively normalized by
the resistance to flow $R_{0}$ and the area $\Lambda_{0}$ in the absence of
gel deformation. (d) A comparison of the measured normalized deformation of
the microgel $l/d$ with its predicted deformation based on solid contact,
$l_{\rm theory}/d$. The graphs include experimental (color) and simulation
results (white).

In the absence of deformation, the pressure drop in the channel occurred both
around the spherical gel particle and in the slit, as shown by the isobars in
the fluid simulations of Fig. 3b (left panel). When the gel deformed however,
the pressure drop was focused over a short distance upstream of the slit,
corresponding to the region where the flow was focused through a narrow gap,
whose cross-sectional area we denote $\Lambda_{\rm exp}$ (Fig. 3b, right
panel). This concentration of the pressure drop indicates that the resistance
to flow was due not to the finger of gel elongating into the slit but rather
to the bulb-like part of the gel that obstructed the fluid flow upstream of
it.

The physical reason behind the focusing of the pressure drop can be understood
by measuring the dependence of $R_{\rm gel}$ on $\Lambda_{\rm exp}$. This was
done by first measuring the area $\Lambda_{\rm sim}$ of this spacing on all
numerical simulations (see inset in Fig. 3c). A geometric argument (see Supp.
Mat. for details) yields a relation between $\Lambda_{\rm sim}$ and other
geometrical parameters:

$\Lambda_{\rm
sim}=\left(\frac{d}{4}-\frac{a}{2}\right)\sqrt{\left(\frac{d}{2}-l-\frac{a}{2}\right)^{2}+\left(\frac{w}{2}+\frac{r_{c}}{\sqrt{2}}\right)^{2}},$
(1)

where $d$ is the bead diameter, $w$ is the width of the slit, $r_{c}$ is the
radius of curvature of the corner of the slit (see Fig. 3a), and $a$ is in
effect the farthest extent of the contact area between the gel and the corner
of the slit (see Fig. 3b, right panel).

Because $a$ is in a direction perpendicular to the imaging plane, it is
inferred by considering the contact mechanics problem of a spherical elastic
bead, of diameter $d$ and elastic modulus $E^{\star}$, being pressed with a
pressure $P_{\rm gel}$ on top of two cylinders of equal radius $r_{c}$ and
spaced by $w$ [35]. The displacement of the gel $l_{\rm theory}$ due to this
forcing was found to relate to the major axis $a$ of the ellipsoidal contact
between the microgel and the trap, modulated by a combination of the
geometrical parameters of the problem (see Supp. Mat. for full derivation):

$l_{\rm theory}=\frac{2a^{2}}{d\sin{\theta}}\mathcal{F}_{1}(r_{c},d),$ (2)

where $\mathcal{F}_{1}(r_{c},d)$ is a function that describes the shape of
contact (see Supp. Mat.) and $\theta$ is the angle created between the
microgel, the corner and the slit (see Fig. 3a).

The experimental area $\Lambda_{\rm exp}$ of the spacing between gel bead and
sidewalls was similarly estimated from Eq. (1), by measuring the geometric
parameters $d,w,r_{c}$ and $l$. Using Eq. (2), we were able to estimate $a$ at
each imposed flow rate as well.

Given the above geometric measurements it is now possible to determine the
scaling of $R_{\rm gel}$ with the interstitial space $\Lambda$. Calling
$R_{0}$ and $\Lambda_{0}$ the values of the parameters in the absence of gel
deformation, _i.e._ at the lowest flow rate tested, we plot $R_{\rm
gel}/R_{0}$ as a function of $(\Lambda_{0}/\Lambda_{\rm exp})^{2}$ in Fig. 3c.

The plot shows an excellent collapse for all experimental and numerical
results and indicates that $R_{\rm gel}$ scales as $\Lambda_{\rm exp}^{-2}$.
The collapse of the measurements of $R_{\rm gel}$ on a single master curve
shows that the resistance to flow is indeed due to the deformation of the soft
solid upstream of the slit which in turn determines the size of the gap that
the flow must go through. Moreover, the rapid increase of $R_{\rm gel}$ makes
it the dominant source of pressure drop compared with a the other sections in
the rest of the microchannel.

In turn the added resistance due to flow focusing couples back to modify the
shape of the microgel. This determines the values of $l$ and $a$, which are
related together by Eq. (2). Contact mechanics modeling shows that $a$ depends
on the ratio $P_{gel}/E^{\star}$ and a combination of the geometric
parameters, which yields an implicit relationship between mechanical and
geometric effects:

$\displaystyle\left(\frac{a}{d}\right)^{3}=$ $\displaystyle\frac{3P_{\rm
gel}}{E^{\star}}\;\times$
$\displaystyle\left[\frac{aw+\frac{2}{3}(\frac{d}{2}-a)[w+r_{c}(1-\cos{\theta})]}{4\pi
d^{2}\sin{\theta}}\right]\mathcal{F}_{2}(r_{c},d)$ (3)

where $\mathcal{F}_{2}(r_{c},d)$ is another function of the shape of contact
with the corner (see Supp. Mat.).

The expected theoretical value of the rescaled gel elongation $l_{\rm
theory}/d$ was calculated from Eqs. (2) and (3) and compared to the
experimentally and numerically computed values of elongation, as shown in Fig.
3d. Again the theory collapses the data onto a master curve, particularly for
small deformations for which contact mechanics are expected to apply.

It is now possible to quantitatively understand the relationship between the
pressure $P_{\rm gel}$ on the trapped soft bead and the flow $Q_{\rm thru}$ it
lets pass through. Combining Eqs. (1)-(3), we find the following relationship
between the normalized flow rate through the slit $\mathcal{Q}={Q_{\rm
thru}R_{0}}/{E^{\star}}$ and the normalized pressure drop $\mathcal{P}=P_{\rm
gel}/E^{\star}$ across the gel:

$\mathcal{Q}=C_{1}\mathcal{P}\left(1-C_{2}\mathcal{P}^{\frac{1}{3}}\right)^{2}\left(\left(1-C_{3}\mathcal{P}^{\frac{2}{3}}\right)^{2}+C_{4}^{2}\right)$
(4)

where the $C_{i}$ are geometric parameters dependent upon $w/d$ and $r_{c}/d$
(see Supp. Mat.).

From Eq. (4) we can predict the evolution of $Q_{\rm thru}$ as the driving
pressure is increased. Experiment and theory are in good agreement, see Fig.
4a. For small values of $\mathcal{P}$ the pressure does not lead to
significant deformation of the microgel, such that further increasing
$\mathcal{P}$ leads to a nearly linear increase in $\mathcal{Q}$. When the
pressure is increased beyond $\mathcal{P}\approx 1$ ($P_{\rm gel}\approx
E^{\star}$) the gel bead deforms, which increases the value of $R_{\rm gel}$
and leads to a decrease in the flow rate through the slit. The transition
between the two regimes depends on the slit geometry and shifts to lower
values of $\mathcal{P}$ for larger values of $w$, indicating that wider slits
are easier to plug than thinner slits.

Figure 4:  (a) Comparison of flow rate through the obstructed thrupass
channel, shown in points, with the predicted flow rate, shown as dashed lines.
Normalized flow rates, $\mathcal{Q}$, are given as a function of the
normalized pressure $\mathcal{P}$ for the full range of pressures tested
experimentally. (b) Maximum flow rate $Q_{\rm max}$ past the microgel as a
function of the Young’s modulus $E^{\star}$ of the trapped bead. Points:
Experimental data. Error bars represent 1 standard deviation. Line shown
represents the linear best fit.

Finally, Eq. (4) shows that the dimensionless flow rate $\mathcal{Q}$ depends
only on the dimensionless pressure $\mathcal{P}$ and on the geometry of the
slit, through the parameters $C_{1}-C_{4}$. For a given geometry it follows
that $\mathcal{Q}$ traces a unique curve, of parameter $\mathcal{P}$, whose
maximum value $\mathcal{Q}_{\rm max}$ thus only depends on the geometry of the
slit. Therefore the maximum dimensional flow rate ${Q}_{\rm
max}=\mathcal{Q}_{\rm max}E^{\star}/R_{0}$ scales linearly with the Young
modulus of the gel $E^{\star}$: a soft gel deforms right away and plugs the
channel at low pressures, while a stiff gel allows higher flow to go through.
This linear increase is confirmed by comparing the largest measured flow rate
with the prediction of $Q_{\rm max}$, as shown in Fig. 4b.

The results presented here can impact several applications, beyond their
scientific interest. They provide a physical basis to understand the
encapsulation of soft beads in droplets, which has emerged as an important
microfluidic technology [36]. In a different operation regime, the strong
nonlinear relationship between pressure and flow rate can lead to the design
of microfluidic nonlinear flow elements, such as check-valves or flow
limiters. These devices play an important role in ensuring the robustness of
fluidic circuits and protecting against surges. The analysis above shows that
the maximum allowable flow rate scales linearly with $E^{\star}$, thus
providing a simple design rule. Finally, the design of the slit presented here
can also serve to measure the elastic modulus of soft materials, similarly to
a micropipette aspiration device.

## Acknowledgements:

The authors acknowledge microfabrication assistance of Caroline Frot at Ecole
Polytechnique and the Biomaterials and Microfluidics platform platform at
Institut Pasteur. Useful discussions with Hiba Belkadi and Mahdi Daei Daei are
also acknowledged.

## References

  * Païdoussis  _et al._ [2010] M. P. Païdoussis, S. J. Price, and E. De Langre,  _Fluid-structure interactions: cross-flow-induced instabilities_ (Cambridge University Press, 2010).
  * Duprat and Shore [2015] C. Duprat and H. A. Shore,  _Fluid-structure interactions in low-Reynolds-number flows_ (Royal Society of Chemistry, 2015).
  * Liu  _et al._ [2018] Y. Liu, B. Chakrabarti, D. Saintillan, A. Lindner, and O. Du Roure, Proceedings of the National Academy of Sciences 115, 9438 (2018).
  * Duprat  _et al._ [2015] C. Duprat, H. Berthet, J. S. Wexler, O. Du Roure, and A. Lindner, Lab on a Chip 15, 244 (2015).
  * Du Roure  _et al._ [2019] O. Du Roure, A. Lindner, E. N. Nazockdast, and M. J. Shelley, Annual Review of Fluid Mechanics 51, 539 (2019).
  * Baroud  _et al._ [2006] C. N. Baroud, S. Tsikata, and M. Heil, Journal of Fluid Mechanics 546, 285 (2006).
  * Pihler-Puzović  _et al._ [2013] D. Pihler-Puzović, R. Périllat, M. Russell, A. Juel, and M. Heil, Journal of Fluid Mechanics 731, 162 (2013).
  * Hazel and Heil [2005] A. L. Hazel and M. Heil, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 461, 1847 (2005).
  * Jensen  _et al._ [2016] K. H. Jensen, K. Berg-Sørensen, H. Bruus, N. M. Holbrook, J. Liesche, A. Schulz, M. A. Zwieniecki, and T. Bohr, Reviews of Modern Physics 88, 035007 (2016).
  * Kim and Beebe [2007] D. Kim and D. J. Beebe, Sensors & Actuators: A. Physical 136, 426 (2007).
  * Louf  _et al._ [2020] J.-F. Louf, J. Knoblauch, and K. H. Jensen, Physical review letters 125, 098101 (2020).
  * Chappel [2020] E. Chappel, Applied Sciences 10, 8858 (2020).
  * Abate  _et al._ [2009] A. R. Abate, C.-H. Chen, J. J. Agresti, and D. A. Weitz, Lab on a chip 9, 2628 (2009).
  * Klein  _et al._ [2015] A. M. Klein, L. Mazutis, I. Akartuna, N. Tallapragada, A. Veres, V. Li, L. Peshkin, D. A. Weitz, and M. W. Kirschner, Cell 161, 1187 (2015).
  * Luo  _et al._ [2014] Y. Luo, D. Chen, Y. Zhao, C. Wei, X. Zhao, W. Yue, R. Long, J. Wang, and J. Chen, Sensors and Actuators B: Chemical 202, 1183 (2014).
  * Elias  _et al._ [2020] M. Elias, A. Dutoya, A. Laborde, A. Lecestre, C. Montis, L. Caselli, D. Berti, B. Lonetti, C. Roux, and P. Joseph, Micro and Nano Engineering 8, 100064 (2020).
  * Preira  _et al._ [2013] P. Preira, M.-P. Valignat, J. Bico, and O. Théodoly, Biomicrofluidics 7, 024111 (2013).
  * Tlili  _et al._ [2022] S. Tlili, F. Graner, and H. Delanoë-Ayari,  _A microfluidic platform to investigate the role of mechanical constraints on tissue reorganization_, preprint (Biophysics, 2022).
  * Zhang  _et al._ [2018] Z. Zhang, J. Xu, and C. Drapaca, Microfluidics and Nanofluidics 22, 120 (2018).
  * Xu  _et al._ [2022] Y. Xu, H. Zhu, Y. Shen, A. P. M. Guttenplan, K. L. Saar, Y. Lu, D. Vigolo, L. S. Itzhaki, and T. P. J. Knowles, MRS Bulletin 47, 119 (2022).
  * Harth  _et al._ [2020] K. Harth, J. Wang, T. Börzsönyi, and R. Stannarius, Soft Matter 16, 8013 (2020).
  * Dressaire and Sauret [2017] E. Dressaire and A. Sauret, Soft Matter 13, 37 (2017).
  * Agbangla  _et al._ [2014] G. C. Agbangla, P. Bacchin, and E. Climent, Soft Matter 10, 6303 (2014).
  * Tao  _et al._ [2021] R. Tao, M. Wilson, and E. R. Weeks, Phys. Rev. E 104, 044909 (2021).
  * Hong  _et al._ [2022] X. Hong, K. W. Desmond, D. Chen, and E. R. Weeks, Phys. Rev. E 105, 014603 (2022).
  * Alborzi  _et al._ [2022] S. Alborzi, B. Clark, and S. M. Hashmi, Soft Matter (2022).
  * Yoon  _et al._ [2016] Y. Yoon, S. Kim, J. Lee, J. Choi, R.-K. Kim, S.-J. Lee, O. Sul, and S.-B. Lee, Scientific Reports 6, 26531 (2016).
  * Khan  _et al._ [2017] Z. S. Khan, N. Kamyabi, F. Hussain, and S. A. Vanapalli, Convergent Science Physical Oncology 3, 024001 (2017).
  * Thielicke and Stamhuis [2014] W. Thielicke and E. J. Stamhuis, Journal of Open Research Software 2, doi:10.5334/jors.bl. (2014).
  * Guermonprez  _et al._ [2015] C. Guermonprez, S. Michelin, and C. N. Baroud, Biomicrofluidics 9, 054119 (2015).
  * Vanapalli  _et al._ [2007] S. A. Vanapalli, D. van den Ende, M. H. G. Duits, and F. Mugele, Applied Physics Letters 90, 114109 (2007).
  * Dendukuri  _et al._ [2005] D. Dendukuri, K. Tsoi, T. A. Hatton, and P. S. Doyle, Langmuir 21, 2113 (2005).
  * Guillou  _et al._ [2016] L. Guillou, A. Babataheri, P.-H. Puech, A. I. Barakat, and J. Husson, Scientific Reports 6, 21529 (2016).
  * Bruus [2007] H. Bruus,  _Theoretical microfluidics_ , Vol. 18 (Oxford university press, 2007).
  * Barber [2018] J. Barber,  _Contact Mechanics_, Solid Mechanics and Its Applications, Vol. 250 (Springer International Publishing, Cham, 2018).
  * Zilionis  _et al._ [2016] R. Zilionis, J. Nainys, A. Veres, V. Savova, D. Zemmour, A. M. Klein, and L. Mazutis, Nature protocols 12, 44 (2016).