# Weakening the effect of boundaries: ‘diffusion-free’ boundary conditions as
a ‘do least harm’ alternative to Neumann

Yufeng Lin${{\dagger}}$ and Rich R. Kerswell${{\ddagger}}$ ∗
${{\dagger}}$ Department of Earth and Space Sciences ∗ Corresponding author.
Email<EMAIL_ADDRESS>Southern University of Science and Technology  Shenzhen
518055  P.R. China
${{\ddagger}}$DAMTP  Centre for Mathematical Sciences  University of Cambridge
CB3 0WA  UK

###### Abstract

In this note, we discuss a poorly known alternative boundary condition to the
usual Neumann or ‘stress-free’ boundary condition typically used to weaken
boundary layers when diffusion is present but very small. These ‘diffusion-
free’ boundary conditions were first developed (as far as the authors know) in
1995 (Sureshkumar & Beris J. Non-Newtonian Fluid Mech. 60, 53-80, 1995) in
viscoelastic flow modelling but are worthy of general consideration in other
research areas. To illustrate their use, we solve two simple ODE problems and
then treat a PDE problem - the inertial wave eigenvalue problem in a rotating
cylinder, a sphere and spherical shell for small but non-zero Ekman number
$E$. Where inviscid inertial waves exist (cylinder and sphere), the viscous
flows in the Ekman boundary layer are $O(E^{1/2})$ weaker than for the
corresponding stress-free layer and fully $O(E)$ weaker than in a non-slip
layer. This reduction could allow precious numerical resources to focus on
other areas of the flow and thereby make smaller, more realistic values of
diffusion accessible to simulations.

###### keywords:

##  1 Introduction

In many flow situations of interest, diffusion of a given physical quantity
may be so small relative to other effects that it would appear safely
neglectable. The well-known caveat to this is the presence of boundaries where
small scales can make diffusion a leading order effect to accommodate the
boundary conditions. The situation can be even more complicated if these
boundary layers also generate internal shear layers in which diffusion plays
an enhanced role in the bulk (e.g. the inertial wave problem to be discussed
below). The strategy almost always is either to take the smallest diffusion
coefficient possible which still allows accurate solutions to be computed
(even though this may be orders of magnitude bigger than the actual value)
and/or the boundary conditions are adjusted to weaken the boundary layers and
possibly internal layers which result. In the case of high Reynolds number
flows (or low Ekman numbers) where momentum diffusion is small, the latter is
invariably achieved by adopting ‘stress-free’ boundary conditions instead of
the usual ‘non-slip’ conditions of the flow at a solid boundary (e.g. Vidal &
Cebron, 2023). This, however, can create its own problems in that stress-free
boundary conditions prevent any angular momentum transfer through viscous
stresses in rotating systems (Vidal & Cebron, 2023).

The purpose of this short note is to describe a further option - hereafter
referred to as ‘diffusion-free’ boundary conditions - for mitigating the
effects of diffusion and to demonstrate that they produce an even weaker
viscous boundary layer correction than ‘stress-free’ or Neumann conditions. In
the interests of clarity, we specialise below to the vanishing viscosity
situation although it should be clear that what follows is perfectly general
for handling the diffusion of any physical quantity, e.g. a scalar like
density, another vector such as a magnetic field or even a tensor like the
conformation tensor in viscoelastic fluids. In fact, it is within the
viscoelastic flow community that the new boundary condition was (as far as the
authors are aware) first used (Sureshkumar & Beris, 1995; Sureshkumar et al.,
1997). Here, the correct boundary condition on the conformation tensor, which
describes an ensemble-averaged configuration of the polymers, is unclear and
so those authors decided to simply turn the diffusion off at the boundary.
This meant that the diffusion-free equation was imposed at the wall as the
required extra boundary conditions or, equivalently, for a $C^{2}$ solution,
the (tensorial) Laplacian of the conformation tensor had to vanish. Carrying
this over to the flow velocity, the idea is to augment the usual (inviscid)
no-penetration normal condition with the extra conditions that the tangential
components of the inviscid equation hold there. Again for a sufficiently
smooth velocity field, this is equivalent to imposing that the tangential
components of the viscous (vectorial) Laplacian term vanish. Imposing that the
second derivative must vanish when solving a second order equation appears
paradoxical. But, in fact, all this does for a smooth solution is impose that
the ‘rest’ of the equation (comprising of first and zeroth derivative terms)
vanishes at the boundary so the diffusion-free boundary condition is
equivalent to a generalised Robin condition (‘generalised’ as this condition
can be a nonlinear relation between zeroth and first derivatives).

Boundary conditions which include the second derivative for a second order
equation are not new. They have been been studied since the 1950s (Feller,
1952; Ventcel, 1959) in the mathematical analysis literature and variously go
under the names of ‘Ventcel’ (e.g. Buffe (2017)), ‘Ventsell’ (e.g. Hinz et al.
(2018)) or ‘Wentzell’ (e.g. Denk et al. (2021)) boundary conditions. The
novelty introduced by Sureshkumar & Beris (1995) was to implicitly use a
special form of them when artificially smoothening a differential equation.

We first illustrate the diffusion-free boundary condition for a couple of
simple ODE problems and then to demonstrate how it works for a realistic PDE
case. The inertial wave problem is chosen for the latter as: a) it is relevant
for rapidly rotating systems such as the Earth’s core (Vidal & Cebron, 2023)
where the exact tangential boundary conditions are unclear; b) it has some
challenging viscous features to resolve in the form of internal shear layers
(McEwan, 1970); c) it is still an area of active research (e.g. He et al.,
2022, 2023; Le Dizes, 2023); and finally d) it is an opportunity to update
some 30-year-old work by one of the authors (Kerswell & Barenghi, 1995,
hereafter KB95).

##  2 Model Problem 1: Nonlinear Robin condition

We first start by considering the simple, albeit nonlinear, ODE problem

$y^{\prime}=-y^{2}$ (1)

over $x\in[0,1]$ ($y^{\prime}:=dy/dx$) with boundary condition $y(0)=1$. This
has the solution $y=1/(1+x)$. Now imagine that there is actually very weak
diffusion in the system and really the equation is

$y^{\prime}=-y^{2}+\varepsilon y{{}^{\prime\prime}}$ (2)

with the diffusion coefficient $\varepsilon\ll 1$. The question then arises as
to what extra boundary condition to impose at $x=1$ because the problem is now
second order. If known, the physically-relevant condition is the natural one,
of course, but sometimes this is unclear or there is a need to minimise the
effect of the boundary to examine its influence or perhaps reduce the
numerical resolution needed to capture the interior solution. Then there is
interest in selecting the boundary condition which does ’least harm’ to the
interior solution. To explore this, we now examine the effect of imposing: (i)
a Dirichlet condition $y(1)=0$; (ii) a Neumann condition $y^{\prime}(1)=0$; or
(iii) the new ‘diffusion-free’ condition $y^{\prime}(1)=-y(1)^{2}$ which is
equivalent to $y^{\prime\prime}(1)=0$ for a smooth solution.

The three problems can be straightforwardly solved using matched asymptotic
expansions to give the leading composite solutions:

$\begin{array}[]{lll}(i)&y=\displaystyle{\frac{1}{1+x}-\frac{1}{2}e^{-(1-x)/\varepsilon}}&{\rm
Dirichlet},\\\ &&\\\ (ii)&y=\displaystyle{\frac{1}{1+x}+\frac{1}{4}\varepsilon
e^{-(1-x)/\varepsilon}}&{\rm Neumann},\\\ &&\\\
(iii)&y=\displaystyle{\frac{1}{1+x}-\frac{1}{4}\varepsilon^{2}e^{-(1-x)/\varepsilon}}&{\rm
Diffusion\\!\\!-\\!\\!free}.\end{array}$ (3)

The first term on the RHS represents the leading ‘outer’ solution (with the
next correction being $O(\varepsilon)$) and each second term is the leading
‘inner’ boundary layer correction. The take-home message here is that this is
$O(1)$ for the Dirichlet condition, $O(\varepsilon)$ for the Neumann condition
and $O(\varepsilon^{2})$ for the diffusion-free condition which therefore
creates the weakest boundary layer correction. It’s worth remarking also that
the actual condition imposed is $y^{\prime}(1)=-y(1)^{2}$ which is not the
usual (linear) Robin condition due to the nonlinearity in $y$ but a suitable
generalisation i.e. a nonlinear relationship between $y$ and $y^{\prime}$
which is imposed at $x=1$.

##  3 Model Problem 2: Interior smoothening

The second ODE problem we consider is

$y=(2+x^{2})H(x)\,:=\left\\{\begin{array}[]{lr}2+x^{2}&\qquad x\geq 0\\\
0&\qquad x<0\end{array}\right.$ (4)

for $x\in[-1,1]$ and $H(x)$ is the Heaviside function. This has a
discontinuity at $x=0$ such as could be present in a hyperbolic system and
non-zero first and second derivatives at one boundary for reasons which will
become clear immediately below. To smooth the solution (e.g. to numerically
estimate it), diffusion could be added as follows

$y=(2+x^{2})H(x)+\varepsilon^{2}y^{\prime\prime}$ (5)

to target the discontinuity hopefully leaving the rest of the solution largely
unchanged (the diffusion coefficient $\varepsilon^{2}$ controls the amount of
smoothening). This makes the problem second order requiring two new boundary
conditions. We now show again that the effect of the diffusion-free boundary
conditions will be less disruptive on the original solution than adopting
Neumann (‘stress-free’) conditions.

Figure 1:  Left: The difference between the stress-free solution $y_{s}$ and
the diffusion-free solution $y_{0}$ is plotted on a log scale in red together
with the corresponding difference between the diffusion-free boundary
condition solution (9), $y_{d}$, and $y_{0}$ (in blue) for $\varepsilon=0.01$
near $x=1$. This plot demonstrates that the stress-free solution produces a
larger deviation at the boundary in $y$: both have an unavoidable interior
deviation of $O(\varepsilon^{2})$ due to the change in the ODE and an $O(1)$
deviation at the discontinuity where $x=O(\varepsilon)$. Right: a numerical
solution for (9) shown as the thin blue line again for $\varepsilon=0.01$
using a spectral method with 300 Chebyshev polynomials (blue dots indicate the
collocation points used; it is expected that $O(1/\varepsilon)$ polynomials
are needed due to sparse collocation point spacing at the middle of the
$[-1,1]$ range). The thick transparent red line is the analytic expression
$y_{d}(x)$ given in (9). Note that the same result is obtained imposing the
equivalent boundary condition that $y^{\prime\prime}=0$ at $x=\pm 1$.

###  3.1 Neumann boundary conditions

Neumann (stress-free) boundary conditions are

$y^{\prime}(-1)=0\quad\&\quad y^{\prime}(1)=0$ (6)

which are often viewed as minimally disruptive. This choice leads to the
leading asymptotic solution

$y_{s}\sim\left\\{\begin{array}[]{cr}\displaystyle{2(1+\varepsilon^{2})e^{-1/\varepsilon}\cosh\left(\frac{1+x}{\varepsilon}\right)}&x<0,\\\
\displaystyle{2+x^{2}+2\varepsilon^{2}+2\varepsilon\sinh\left(\frac{1-x}{\varepsilon}\right)-\biggl{[}\,2\varepsilon+2(1+\varepsilon^{2})e^{-1/\varepsilon}\,\biggr{]}\cosh\left(\frac{1-x}{\varepsilon}\right)}&x\geq
0.\end{array}\right.$ (7)

There are then 3 distinct scalings for how this solution differs from the
diffusion-free solution $y_{0}$ of (4): for $x=0(\varepsilon)$,
$y_{s}-y_{0}=O(1)$, for $1-x=O(\varepsilon)$, $y_{s}-y_{0}=O(\varepsilon)$ to
fix up the boundary condition, and otherwise $y_{s}-y_{0}=O(\varepsilon^{2})$
\- see Figure 1.

###  3.2 Diffusion-free boundary conditions

We now impose ‘diffusion-free’ boundary conditions which are just

$y(-1)=(2+x^{2})H(x)|_{x=-1}=0\quad\&\quad y(1)=(2+x^{2})H(x)|_{x=1}=3$ (8)

The leading asymptotic solution is then

$y_{d}\sim\left\\{\begin{array}[]{cr}\displaystyle{2(1+\varepsilon^{2})e^{-1/\varepsilon}\sinh\left(\frac{1+x}{\varepsilon}\right)}&x<0,\\\
\displaystyle{2+x^{2}+2\varepsilon^{2}\biggl{[}\,1-\cosh\left(\frac{1-x}{\varepsilon}\right)\,\biggr{]}+\biggl{[}\,2\varepsilon^{2}-2(1+\varepsilon^{2})e^{-1/\varepsilon}\,\biggr{]}\sinh\left(\frac{1-x}{\varepsilon}\right)}&x\geq
0.\end{array}\right.$ (9)

There are now only 2 distinct scalings for how this solution differs from the
solution $y_{0}$ of (4): for $x=O(\varepsilon)$, $y_{d}-y_{0}=O(1)$ otherwise
$y_{s}-y_{0}=O(\varepsilon^{2})$ \- see Figure 1. The key point here is that
the order of the correction does not change at the boundary in contrast to the
Neumann case where the correction grows to $O(\varepsilon)$. In fact, in this
model problem the correction goes to precisely zero there because the
$\varepsilon=0$ problem has a well defined boundary value. Hence, again,
diffusion-free conditions are better ‘do least harm’ conditions than Neumann.

The problem (5) & (8) can also be solved numerically using a standard
Chebyshev expansion and shows no issues - see Figure 1(right). It’s worth
remarking that the same solution is obtained (up to discretization errors) if
the boundary conditions $y^{\prime\prime}(\pm 1)=0$ are used instead. These
are equivalent for a smooth solution of course.

##  4 Inertial wave problem in a cylinder

We now switch our attention to a more serious PDE example. The linearised
equations for describing small motions of a viscous, incompressible fluid away
from uniform rotation are

$\frac{\partial\mathbf{u}}{\partial
t}+2\mathbf{\hat{z}}\times\mathbf{u}+\mbox{\boldmath$\nabla$}p=E\nabla^{2}\mathbf{u},$
(10) $\mbox{\boldmath$\nabla$}.\mathbf{u}=0,$ (11)

in the rotating frame. This is a core problem in rotating flows with
appropriately more complicated analogues when further physics is present (e.g.
stratification and/or magnetic fields). It also displays a number of
fascinating viscous features typified by boundary layer breakdowns which give
rise to internal shear layers and a host of different viscous lengthscales
(e.g. see figure 4 of He et al. (2022)). The basic rotation rate $\Omega$ and
cylindrical radius $S$ have been used to non-dimensionalise the system, with
the Ekman number $E=\nu/\Omega S^{2}$ appearing as the non-dimensionalization
of the kinematic viscosity $\nu$.

The inviscid limit, in which $E$ is set to zero and only a no-normal velocity
boundary condition needs be applied, gives rise to the well known inertial
wave problem, where wave solutions

$[\mathbf{u},p]=[u(s,z),v(s,z),w(s,z),\Phi(s,z)]e^{i(m\phi+\lambda t)},$

may be sought (e.g. see Greenspan (1968); Zhang & Liao (2017)), and
$(s,\phi,z)$ are cylindrical coordinates). The reduction of the problem to one
for only the pressure realises the Poincaré (1910) equation

$\frac{1}{s}\frac{\partial}{\partial s}s\frac{\partial\Phi}{\partial
s}-\frac{m^{2}}{s^{2}}\Phi+\left(1-\frac{4}{\lambda^{2}}\right)\frac{\partial^{2}\Phi}{\partial
z^{2}}=0$ (12)

to be solved subject to the boundary conditions

$\displaystyle\frac{\partial\Phi}{\partial s}+\frac{2m}{\lambda s}\Phi$
$\displaystyle=$ $\displaystyle 0\quad{\rm on}\,s=1,$
$\displaystyle\frac{\partial\Phi}{\partial z}$ $\displaystyle=$ $\displaystyle
0\quad{\rm on}\,z=0,\,d.$

For the cylinder problem, Kelvin (1880) first wrote down the separable
solutions

$\mathbf{u}_{mnl}=\frac{e^{i(m\phi+\lambda
t)}}{2(4-\lambda^{2})}\left[\begin{array}[]{c}i\\{(\lambda+2)J_{m-1}(ks)-(\lambda-2)J_{m+1}(ks)\\}\cos(l\pi
z/d)\\\ -\\{(\lambda+2)J_{m-1}(ks)+(\lambda-2)J_{m+1}(ks)\\}\cos(l\pi z/d)\\\
2i\lambda kd(\pi l)^{-1}J_{m}(ks)\sin(l\pi z/d)\end{array}\right],$ (13)
$p=-\frac{1}{k}J_{m}(ks)\cos(l\pi z/d)e^{i(m\phi+\lambda t)},$ (14)

where

$\lambda=\frac{\pm 2}{\sqrt{1+k^{2}d^{2}/(\pi l)^{2}}},$

and $k$ is a solution, indexed by $n$ such that
$0\,<\,k_{n=1}\,<\,k_{n=2}\ldots$, of

$s\frac{d}{ds}J_{m}(ks)+\frac{2m}{\lambda}J_{m}(ks)=0|_{s=1}.$

Three numbers, $l$, $n\in{\mathbb{N}}$, and $m\in{\mathbb{Z}}$, are sufficient
to specify the inertial mode and in particular the frequency
$\lambda=\lambda_{lmn}$. These indices correspond roughly with the number of
nodes axially, radially, and azimuthally respectively in the pressure
eigenfunction. As in KB95 and Kerswell (1999), we focus here on the inertial
mode with $m=n=l=1$ and $d=1.9898$ (the quoted value in KB95) with inertial
velocity field $\mathbf{u}_{0}$ and real frequency
$\lambda_{0}=1.000007648337\approx 1$. The subscript is to indicate that these
are the $E=0$ quantities compared to those, $\mathbf{u}_{E}$ and
$\lambda_{E}$, for $E\neq 0$. The viscous correction is then defined as
$\mathbf{u}_{c}:=\mathbf{u}_{E}-\mathbf{u}_{0}$. Reinstating viscosity
requires 2 extra (tangential) boundary conditions to supplement to no-
penetration (normal) condition.

###  4.1 Non-slip boundary conditions

If the container walls are rigid then non-slip boundary conditions are
required which is the usual case discussed in text books (e.g. §2.5 in
Greenspan (1968)). If $\mathbf{u}_{0}(\mathbf{x})e^{i\lambda t}$ is the
solution to the inertial wave problem (12) with $\lambda$ a real frequency in
the interval $[-2,2]$ for $E=0$, then the leading viscous correction for
$E\neq 0$ is the $O(E^{0})$ boundary velocity
${\mathbf{\tilde{u}}}_{0}(\mathbf{x},t)$ forced by the non-slip condition

$\mathbf{u}_{0}+{\mathbf{\tilde{u}}}_{0}={\bf 0}|_{\partial V}.$ (15)

(hereafter $\mathbf{\tilde{\,\,}}$ indicates a boundary layer flow). This
decays into the interior from the boundary over an Ekman boundary of thickness
$E^{1/2}$ but drives an $O(E^{1/2})$ viscous flow
$\mathbf{u}_{1}(\mathbf{x},t)$ in the interior by ‘Ekman pumping’ (equation
2.6.13 in Greenspan (1968)). These flows mean that $\lambda$ becomes complex
with a viscous decay rate $\lambda_{i}:=\Im m(\lambda)=O(E^{1/2})\,>\,0$.
There are also shear layers spawned by the corners of the cylinder (e.g.
McEwan (1970) and figures 2-4 below) but these appear unimportant for the
leading decay rate calculation (§2.9 in Greenspan (1968)).

###  4.2 Stress-free boundary conditions

For other situations, the boundary may be stress-free (e.g. for a gaseous
planet) or this choice is assumed to make the viscous correction at the
boundary as weak as possible (e.g. Lin & Ogilvie (2021); Vidal & Cebron
(2023)). Then the leading viscous correction is
$E^{1/2}{\mathbf{\tilde{u}}}_{1}(\mathbf{x},t)$ in a boundary layer of
thickness $E^{1/2}$ such that, at, for example, the cylinder sidewall,

$s\frac{\partial}{\partial
s}\biggl{(}\frac{v_{0}}{s}\biggr{)}-\frac{\partial\tilde{v}_{1}}{\partial\xi}=0\biggl{|}_{s=1}\biggr{.}\quad\&\quad\frac{\partial
w_{0}}{\partial
s}-\frac{\partial\tilde{w}_{1}}{\partial\xi}=0\biggl{|}_{s=1}\biggr{.}$ (16)

where $\xi:=(1-s)/E^{1/2}$ is the rescaled normal coordinate. The corrective
viscous bulk flow (as opposed to localised shear layers) is then $O(E)$ driven
by the subsequent $O(E)$ Ekman pumping at the boundaries and the bulk viscous
term also of $O(E)$. As a result the viscous decay rate is much smaller -
$\lambda_{i}:=O(E)\,>\,0$.

Unfortunately, for our purposes here, the cylinder has a singular feature for
these boundary conditions: the inertial mode automatically satisfies stress-
free conditions on the top and bottom boundaries, i.e. $\partial
u_{0}/\partial z=\partial v_{0}/\partial z=0$ at $z=0$ and $z=d$ already.
Hence a viscous boundary layer only forms on the side wall and in particular
no internal shear layers are formed by the corner region where the horizontal
boundaries meet the sidewall which is an exceptional response. - e.g. a
calculation in more geophysically-interesting sphere or spherical shell does
have shear layers (e.g. Lin & Ogilvie (2021); Vidal & Cebron (2023) and the
next section). As a result, for stress-free boundary conditions, the
calculation for the viscous decay rate $\lambda_{i}$ in a cylinder is just 1D
- the variation in $z$ can be separated off - and so very small Ekman numbers
can be reached easily (see Table 1).

###  4.3 Diffusion-free boundary conditions

Diffusion-free boundary conditions are

$i\lambda u-2v+\frac{\partial p}{\partial s}=0,\quad i\lambda
v+2u+\frac{imp}{s}=0\quad\&\quad w=0$ (17)

on the top ($z=d$) and bottom ($z=0$) and

$u=0,\quad i\lambda v+2u+\frac{imp}{s}=0\quad\&\quad i\lambda w+\frac{\partial
p}{\partial z}=0$ (18)

on the sides ($s=1$) of the cylinder. For a $C^{2}$ solution this is
equivalent to setting the tangential components of the vectorial Laplacian to
zero. That is, the conditions

$u=0,$
$\mbox{\boldmath$\hat{\phi}$}.\nabla^{2}\mathbf{u}=\frac{\partial^{2}v}{\partial
s^{2}}+\frac{1}{s}\frac{\partial v}{\partial
s}-\frac{m^{2}+1}{s^{2}}v+\frac{\partial^{2}v}{\partial
z^{2}}+\frac{2imu}{s^{2}}=0,$
$\mathbf{\hat{z}}.\nabla^{2}\mathbf{u}=\frac{\partial^{2}w}{\partial
s^{2}}+\frac{1}{s}\frac{\partial w}{\partial
s}-\frac{m^{2}}{s^{2}}w+\frac{\partial^{2}w}{\partial z^{2}}=0$ (19)

imposed at $s=1$ are equivalent to those in (18) having the same complex
frequencies $\lambda(E)$ as eigenvalues (not shown but verified).

###  4.4 Scalings

The inviscid inertial wave satisfies all the diffusion-free boundary
conditions. However the $O(E)$ bulk viscous term drives a primary $O(E)$
interior flow, $E\mathbf{u}_{2}(\mathbf{x},t)$, which does not, leading to an
$O(E)$ error in the diffusion-free boundary conditions (with no boundary-
normal derivatives). These can be fixed up by a very weak $O(E)$ boundary
correction flow, $E{\mathbf{\tilde{u}}}_{2}(\mathbf{x},t)$. The associated
Ekman pumping drives a secondary bulk flow at $O(E^{3/2})$. Hence the
corrective viscous flow is $O(E)$ in both the bulk and the Ekman boundary
layer. In contrast, the secondary flow is $O(E)$ in the bulk and a larger
$O(E^{1/2})$ in the boundary layer for stress-free boundary conditions. For
non-slip conditions, the secondary flow is $O(E^{1/2})$ in the bulk and $O(1)$
in the boundary layer.

###  4.5 Numerical Implementation

We now demonstrate the scaling predictions discussed above using numerical
calculations over a range of small but finite $E$. As in KB95, we exploit the
symmetries of the inertial waves by only solving the eigenvalue problem
(10)-(11) in the upper ‘quadrant’ $(s,z)\in[0,1]\times[d/2,d]$ of an extended
cylinder $(s,z)\in[-1,1]\times[0,d]$. If $m$ (the azimuthal wavenumber) is
even (odd), then $(w,p)$ are even (odd) and $(u,v)$ are odd (even) functions
of $s$. Likewise about the cylindrical equator $z=d/2$, if $l$ is even (odd),
then $(u,v,p)$ are even (odd) and $w$ is odd (even) functions of
$z^{*}:=2z/d-1$. Hence, for the situation of interest here, $m=l=1$, the
appropriate spectral expansions used are

$(u,v)(s,z)=\sum^{N}_{i=1}\sum_{j=1}^{N}(u_{ij},v_{ij})T_{2i-2}(s)T_{2j-1}(z^{*}),$
(20)
$w(s,z)=\sum^{N}_{i=1}\sum_{j=1}^{N}w_{ij}T_{2i-1}(s)T_{2j-2}(z^{*}),\qquad
p(s,z)=\sum^{N}_{i=0}\sum_{j=0}^{N}p_{ij}T_{2i-1}(s)T_{2j-1}(z^{*})$ (21)

where $T_{n}(x):=\cos(n\cos^{-1}x)$ is the $n^{th}$ Chebyshev polynomial. The
equations were collocated over the positive zeros of $T_{2N}$ in $s$ and
$z^{*}$ (the expansion lengths in $r$ and $z$ are taken equal for simplicity).
Inverse iteration was then used to find the relevant (complex) eigenvalue and
eigenfunction as the inviscid real frequency $\lambda_{0}$ provides a good
guess. KB95 talk about a maximal (given the 300 MBytes available) $N=35$ full
eigenvalue solve taking 3 days on a Sun SparcCenter 2000 (the University of
Newcastle’s largest computer in 1995). The equivalent today takes 9.4 mins on
an Apple M1 laptop even using Matlab instead of the Fortran 77 of KB95.
Employing inverse iteration to just focus on the eigenvalue of interest brings
that down to just 18 seconds.

Table 1:  Decay rates $\lambda_{i}$ at various $E$ and different boundary conditions. Non-slip and diffusion-free boundary conditions use up to $N=110$ to check convergence. The stress-free case reduces to a 1D problem and resolution can go much higher: $N=6,000$ and $10,000$ give $\lambda=1.0000076483366+i9.469658E$ for $E=10^{-10}$. Removing the top and bottom boundary conditions for the diffusion-free case gives the starred decay rates on the far right which shows convergence to the interior dissipation rate value of $9.97105E$ indicating that the boundary layer contribution is higher order. $E$ |  | KB95 |  | non-slip |  | stress-free |  | diffusion-free
---|---|---|---|---|---|---|---|---
|  | $\lambda_{i}/\sqrt{E}$ |  | $\lambda_{i}/\sqrt{E}$ |  | $\lambda_{i}/E$ |  | $\lambda_{i}/E$
$10^{-4}$ |  | $1.5665$ |  | $1.5665$ |  | $9.4362$ |  |  $9.831$
$10^{-5}$ |  | $1.4834$ |  | $1.4859$ |  | $9.4591$ |  |  $9.927$
$10^{-6}$ |  |  |  | $1.4619$ |  | $9.4663$ |  |  $9.957$ $9.9647*$
$3\times 10^{-7}$ |  |  |  | $1.4571$ |  |  |  | $9.963$
$10^{-7}$ |  |  |  |  |  | $9.4686$ |  |  $9.9689^{*}$
$10^{-8}$ |  |  |  |  |  | $9.4694$ |  |  $9.9703^{*}$
$10^{-9}$ |  |  |  |  |  | $9.4696$ |  |  $9.9707^{*}$
$10^{-10}$ |  |  |  |  |  | $9.4697$ |  |  $9.9708^{*}$
$10^{-11}$ |  |  |  |  |  |  |  |  $9.9709^{*}$
|  | $\underbrace{\hskip 92.47145pt}_{E\,\rightarrow\,0}$ |  | $\underbrace{\hskip 135.15059pt}_{\text{interior dissipation only}}$
Wedemeyer |  | $1.732$ |  | $9.97105$
Kudlick |  | $1.451$ |  |  |  |
Table 2:  Convergence of decay rate $\lambda_{i}/\sqrt{E}$ estimates with varying $N$ for non-slip boundary conditions ($N=110$ is the largest resolution that runs on a 192GB memory machine). $N$ |  | $E=10^{-4}$ |  | $E=10^{-5}$ |  | $E=10^{-6}$ |  | $E=3\times 10^{-7}$
---|---|---|---|---|---|---|---|---
30 |  | $1.5664732$ |  | $1.4861012$ |  |  |  |
40 |  | $1.5665071$ |  | $1.4870688$ |  |  |  |
50 |  | $1.5665072$ |  | $1.4866554$ |  | $1.4396532$ |  |
60 |  | $1.5665072$ |  | $1.4865893$ |  | $1.4684427$ |  | $1.3328341$
70 |  | $1.5665072$ |  | $1.4865937$ |  | $1.4623715$ |  | $1.4502321$
80 |  |  |  | $1.4865936$ |  | $1.4622908$ |  | $1.4639174$
90 |  |  |  | $1.4865936$ |  | $1.4619161$ |  | $1.4581648$
100 |  |  |  |  |  | $1.4618646$ |  | $1.4571169$
110 |  |  |  |  |  | $1.4618879$ |  | $1.4571179$

Figure 2:  Plots of the viscous corrective field for different boundary
conditions at $E=10^{-4}$: non-slip (left) $\log_{10}|\mathbf{u}_{c}|$ plotted
with $\max|\mathbf{u}_{c}|=4.24$, stress-fress (middle)
$\log_{10}|\mathbf{u}_{c}|/\sqrt{E}$ plotted with
$\max|\mathbf{u}_{c}|=3.03\sqrt{E}$ and diffusion-free conditions (right)
$\log_{10}|\mathbf{u}_{c}|/E$ with $\max|\mathbf{u}_{c}|=44.56E$. This and
figures 3 and 4 show how the viscous corrective flows scale with $E$ for the
various boundary conditions.

###  4.6 Decay rate $\lambda_{i}$ and Eigenfunctions

To identify the viscous corrective flow to the inviscid inertial wave
$\mathbf{u}_{111}$ \- see (13) - due to $E\neq 0$, we use the well-known
orthogonality property of inertial waves (§2.7 Greenspan (1968)) to
renormalise the viscous eigenfunction to have exactly the inertial wave field
$\mathbf{u}_{111}$ and then subtract it -

$\mathbf{u}_{c}:=\biggl{[}\frac{\int_{V}\mathbf{u}^{*}_{111}.\mathbf{u}_{111}\,{\bf
d}^{3}\mathbf{x}}{\int_{V}\mathbf{u}^{*}_{111}.\mathbf{u}_{E}\,{\bf
d}^{3}\mathbf{x}}\biggr{]}\mathbf{u}_{E}-\mathbf{u}_{111}.$ (22)

This definition means that $\mathbf{u}_{c}$ is orthogonal to
$\mathbf{u}_{111}$ and is normalised to be the viscous correction to the field
$\mathbf{u}_{111}$ as defined in (13)

Table 1 collects all the decay rate results across the 3 different boundary
conditions mentioned above. Only the non-slip case was considered in KB95 and
the results there are confirmed for $E=10^{-4}$ and now properly converged for
$10^{-5}$. Treating even smaller $E$ suggests that Kudlick’s formula for the
decay rate is to be preferred to Wedemeyer’s: see KB95 for context. Table 2
shows how the non-slip results converge. The $N=30$ result for $E=10^{-5}$ is
misleadingly good compared to KB95. In fact $N=35$ gives $\lambda_{i}=1.4913$
before convergence really starts at $N=40$.

Figure 3:  Plots of the viscous corrective field for different boundary
conditions at $E=10^{-5}$: non-slip (left) $\log_{10}|\mathbf{u}_{c}|$ plotted
with $\max|\mathbf{u}_{c}|=4.24$, stress-fress (middle)
$\log_{10}|\mathbf{u}_{c}|/\sqrt{E}$ plotted with
$\max|\mathbf{u}_{c}|=3.00\sqrt{E}$ and diffusion-free conditions (right)
$\log_{10}|\mathbf{u}_{c}|/E$ with $\max|\mathbf{u}_{c}|=43.47E$. This and
figures 2 and 4 show how the viscous corrective flows scale with $E$ for the
various boundary conditions.

Figure 4:  Plots of the viscous corrective field for different boundary
conditions at $E=10^{-6}$: non-slip (left) $\log_{10}|\mathbf{u}_{c}|$ plotted
with $\max|\mathbf{u}_{c}|=4.24$, stress-fress (middle)
$\log_{10}|\mathbf{u}_{c}|/\sqrt{E}$ plotted with
$\max|\mathbf{u}_{c}|=2.99\sqrt{E}$ and diffusion-free conditions (right)
$\log_{10}|\mathbf{u}_{c}|/E$ with $\max|\mathbf{u}_{c}|=42.35E$. This and
figures 2 and 3 show how the viscous corrective flows scale with $E$ for the
various boundary conditions.

The contribution to the decay rate from the interior can be estimated as
$(k^{2}+\alpha^{2})E$ (Kerswell (1999)) and with
$\lambda^{2}=4/(1+k^{2}/\alpha^{2})$ is just
$4\alpha^{2}E/\lambda^{2}=9.97104976E$. The diffusion-free decay rate is well
approximated by

$\lambda_{i}=9.97105E-14E^{3/2}+O(E^{2})$ (23)

indicating that the boundary layer contribution is $O(E^{1/2})$ smaller than
the interior’s (this is clearer for the 1D case where the top and bottom
boundaries are made stress-free - see the starred entries in Table 1 which can
be computed to much smaller $E$). This is in contrast to the stress-free case
where the decay rate does not converge to the interior’s contribution as
$E\rightarrow 0$ \- see Table 1 - indicating that the boundary layer
contribution is at the same order as the interior’s. Hence the diffusion-free
boundary condition again is less disruptive than the stress-free (Neumann)
condition.

The various corrective viscous velocities $\mathbf{u}_{c}(r,z)$ are shown in
Figures 2-4 for $E=10^{-4}$, $10^{-5}$ and $10^{-6}$ for the different choices
of boundary conditions. Due to the vastly different amplitudes of these flows,
the leading viscous scaling with $E$ is removed before contouring. That is,
for non-slip boundary conditions, we expect $\mathbf{u}_{c}$ to be $O(1)$ in
the boundary layer so there is no rescaling in this case whereas for the
stress-free case, $\mathbf{u}_{c}$ is $O(E^{1/2})$, so
$\log_{10}|\mathbf{u}_{c}/E^{1/2}|$ is plotted and
$\log_{10}|\mathbf{u}_{c}/E|$ for the diffusion-free case. The various maximum
amplitudes are collected in Table 3 confirming the expected scaling. The
formation of internal shear layers is already evident at $E=10^{-4}$ for the
non-slip and diffusion-free conditions and these get more intense and narrower
as $E$ decreases (their thickness is $E^{1/3}$). Unfortunately for our
purposes here, they are completely absent from the stress-free situation for
reasons already mentioned.

Finally, figure 5 shows the unscaled amplitude of $\mathbf{u}_{c}$ as a
function of $s$ at a given $z$. The top three curves are for the non-slip
situation for $E=10^{-4}$, $10^{-5}$ and $10^{-6}$ and show how the internal
shear layers only really become prominent at $E=10^{-6}$. All have the same
$O(1)$ value at the side wall ($=-|\mathbf{u}_{111}(1,z)|$). The thin black
line is the stress-free profile at $E=10^{-6}$ which has an $O(E^{1/2})$
boundary value and is largely featureless in the interior as already commented
on. Finally the lowest thick red line (computed with $N=90$) and the thin red
line through its centre (computed with $N=110$) are for the diffusion-free
conditions. This is an $O(E^{1/2})$ smaller than the non-slip flow everywhere.
Perhaps, not surprisingly, the velocity profiles are largely similar in shape
as they are formed by the same geometrical effects (the reflections of the
internal shear layers) but the point is, of course, their magnitude is very
different.

Table 3:  Maximum amplitude of the viscous correction velocity $|\mathbf{u}_{c}(r,z)|$ over the cylindrical domain.  $E$ |  | non-slip |  | stress-free |  | diffusion-free
---|---|---|---|---|---|---
$10^{-4}$ |  | $4.24$ |  | $3.03\sqrt{E}$ |  | $44.56E$
$10^{-5}$ |  | $4.24$ |  | $3.00\sqrt{E}$ |  | $43.48E$
$10^{-6}$ |  | $4.24$ |  | $2.99\sqrt{E}$ |  | $42.35E$

Figure 5:  The magnitude of the unscaled corrective viscous velocity,
$|\mathbf{u}_{c}|(s,z^{*})$, plotted as a function of $s$ for a fixed choice
of $z^{*}$ and various different boundary conditions. The upper 3 (blue,
purple and cyan) lines are for non-slip conditions with $z^{*}=d/\sqrt{2}$
(chosen to be arbitrary). The wide transparent blue line corresponds to $N=80$
and the blue line $N=90$ both at $E=10^{-6}$ confirming that the resolution is
sufficient; the purple line is for $E=10^{-5}$ and the cyan line is for
$E=10^{-4}$ (the corresponding dots at $s=1$ are all at $1.32$). The middle
black line is for stress-free conditions at $E=10^{-6}$ ($N=100$) and uses
$z^{*}=d/2$ as the interior happens to rather featureless at
$z^{*}=d/\sqrt{2}\approx 1.4$ \- e.g. see figure 4(middle) (the black dot is
at $2.99\times 10^{-3}$). The lower red curves are for diffusion-free
conditions with $z^{*}=d/\sqrt{2}$ again. The wide transparent red line uses
$N=90$ and the thin red line threading it has $N=110$ (the red dot is at
$1.32\times 10^{-5}$). This plot emphasizes the difference in scales of the
boundary-induced flows between the different boundary conditions and how the
interior structure starts to emerge at $E=10^{-6}$ \- note the structure for
the diffusion-free conditions is very similar to that for the non-slip
conditions albeit 4 orders or magnitude smaller.

##  5 Inertial wave problem in a sphere

The inertial wave problem in a sphere (equations (10) and (11)) also has
separable solutions in the inviscid limit where only the normal velocity needs
to vanish at the walls. These separable solutions involve modified spheroidal
coordinates ($\eta,\mu$) introduced by Bryan (1889) where

$s=\left(\frac{4}{4-\lambda^{2}}-\eta^{2}\right)^{1/2}(1-\mu^{2})^{1/2},\quad
z=\left(\frac{4}{\lambda^{2}}-1\right)^{1/2}\eta\mu.$ (24)

Using these, the pressure field can be written as (e.g. see Greenspan (1968);
Zhang & Liao (2017))

$\Phi_{mnk}=P_{n}^{m}(\eta/c_{mnk})P_{n}^{m}(\mu)\mathrm{e}^{\mathrm{i}(m\phi+\lambda
t)}$ (25)

where $P_{n}^{m}$ represents an associated Legendre function and
$c_{mnk}=(1-\lambda_{mnk}^{2}/4)^{-1/2}$. The real frequency $\lambda_{mnk}$
is the $k$th eigenvalue solution of the transcendental equation

$mP_{n}^{m}(\lambda/2)=2\left(1-\frac{\lambda^{2}}{4}\right)\frac{d}{d\lambda}P^{m}_{n}(\lambda/2),$
(26)

and so each inertial mode is uniquely specified by a triplet of indices
$(m,n,k)$.

###  5.1 Numerical Implementation

The viscous problem ($E\neq 0$) is numerically solved using a pseudo-spectral
method based on an expansion of spherical harmonics on spherical surfaces and
Chebyshev collocation in the radial direction. The incompressible condition
$\nabla\cdot\mathbf{u}=0$ is fulfilled by the the toroidal-poloidal
decomposition

$\mathbf{u}=\mbox{\boldmath$\nabla$}\times(\mathcal{T}\mathbf{\hat{r}})+\mbox{\boldmath$\nabla$}\times\mbox{\boldmath$\nabla$}\times(\mathcal{P}\mathbf{\hat{r}})$
(27)

and, for an inertial wave with the azimuthal wavenumber $m$, the scalars
$\mathcal{T}$ and $\mathcal{P}$ are expanded in spherical harmonics

$\mathcal{T}=\sum_{n=m}^{N}\mathcal{T}_{n}(r)Y_{n}^{m}(\theta,\phi),\qquad\mathcal{P}=\sum_{n=m}^{N}\mathcal{P}_{n}(r)Y_{n}^{m}(\theta,\phi).$
(28)

where $(r,\theta,\phi)$ are spherical polar coordinates. The projections
$\mathbf{\hat{r}}\cdot\mbox{\boldmath$\nabla$}\times$ and
$\mathbf{\hat{r}}\cdot\mbox{\boldmath$\nabla$}\times\mbox{\boldmath$\nabla$}\times$
of the equation

$\mathrm{i}\lambda\mathbf{u}+2\mathbf{\hat{z}}\times\mathbf{u}+\mbox{\boldmath$\nabla$}p=E\nabla^{2}\mathbf{u}$
(29)

onto spherical harmonics, give a set of ODEs for $\mathcal{T}_{n}(r)$ and
$\mathcal{P}_{n}(r)$ with a block tri-diagonal matrix structure. For the
radial dependence, Chebyshev collocation is used on the $K+1$ Gauss-Lobatto
nodes with the projected ODEs replaced by corresponding boundary conditions on
the boundary nodes.

###  5.2 Boundary Conditions

The non-slip boundary conditions $\mathbf{u}={\bf 0}$ correspond to

$\mathcal{T}_{n}(r)=0,\quad\mathcal{P}_{n}(r)=0\quad\&\quad\frac{d\mathcal{P}_{n}(r)}{dr}=0$
(30)

at the boundary $r=1$ (the sphere radius $R$ is used to define the Ekman
number as $E:=\nu/\Omega R^{2}$).

The stress-free boundary conditions are

$u_{r}=0,\quad\frac{\partial}{\partial
r}\left(\frac{u_{\theta}}{r}\right)=0,\quad\frac{\partial}{\partial
r}\left(\frac{u_{\phi}}{r}\right)=0,$ (31)

which translates into

$\mathcal{P}_{n}(r)=0,\quad
2\frac{d\mathcal{P}_{n}(r)}{dr}-r\frac{d^{2}\mathcal{P}_{n}(r)}{dr^{2}}=0,\quad
2\mathcal{T}_{n}(r)-r\frac{d\mathcal{T}_{n}(r)}{dr}=0$ (32)

for the toroidal and poloidal components on $r=1$.

The diffusion-free boundary conditions are

$\mathbf{u}\cdot\mathbf{\hat{r}}=0\quad\&\quad(\mbox{\boldmath$\nabla$}^{2}\mathbf{u})\times\mathbf{\hat{r}}={\bf
0}$ (33)

at $r=1$ or

$\mathcal{P}_{n}(r)=0,\quad
r^{2}\frac{d^{3}\mathcal{P}_{n}(r)}{dr^{3}}-n(n+1)\frac{d\mathcal{P}_{n}(r)}{dr}=0,\quad
r^{2}\frac{d^{2}\mathcal{T}_{n}(r)}{dr^{2}}-n(n+1)\mathcal{T}_{n}(r)=0.$ (34)

At the origin $r=0$, we apply the regularity conditions (Reese et al., 2006)

$\mathcal{P}_{n}=0,\quad\mathcal{T}_{n}=0,$ (35)

for the even degree of spherical harmonics and

$\frac{d\mathcal{P}_{n}}{dr}=0,\quad\frac{d\mathcal{T}_{n}}{dr}=0,$ (36)

for the odd degree of spherical harmonics. The generalized eigenvalue problem
is solved using an iterative method to find the least damped eigenmode around
the inviscid eigenfrequency.

Figure 6:  Decay rate of inertial mode (2,4,1) with the inviscid eigen-
frequency $\lambda_{241}=-1.0926$ in a sphere.

Figure 7:  Viscous corrective field of the full sphere inertial mode (2,4,1)
in the meridional plane for no-slip $|{\bf u}_{c}|$(left), stress-free $|{\bf
u}_{c}|/\sqrt{E}$ (middle), and diffusion-free $|{\bf u}_{c}|/E$ (right)
boundary conditions at $E=10^{-6}$.

Figure 8:  The magnitude of the unscaled corrective viscous velocity along the
dashed lines in Figure 7 for a sphere at $E=10^{-6}$.

###  5.3 Decay rate $\lambda_{i}$ for a sphere

Figure 6 shows the decay rate of a typical inertial mode (here
$(m,n,k)=(2,4,1)$) as a function of the Ekman number for the different
boundary conditions. We find $O(E^{1/2})$ decay rates for non-slip; $O(E)$ for
stress-free and $O(E^{3/2})$ for diffusion-free b.c.s. This indicates that the
viscous correction in the diffusion-free situation is weakest.

Since there are $O(E)$ corrective flows driven in the interior, this last
result looks anomalous until one realises that the interior dissipation
contribution vanishes identically for a sphere (Zhang et al., 2001) and a
spheroid (Zhang et al., 2004). This is easily understood after realising that
the inertial waves in these geometries are polynomials in (the cylindrical
radius) $s$ and $z$ (see the discussion on p118 of Kerswell (1993)) and
partitioned into orthogonal vector spaces ${\mathscr{V}}_{n}$ where $n$ is the
order of the polynomial (see p125 in Kerswell (1993)). The interior
dissipation integral is then

$\int\mathbf{u}_{0}^{*}\cdot\nabla^{2}\mathbf{u}_{0}\,dV=\frac{4}{\lambda^{2}}\int\mathbf{u}_{0}^{*}\cdot\frac{\partial^{2}\mathbf{u}_{0}}{\partial
z^{2}}\,dV=0$ (37)

as $\partial^{2}\mathbf{u}_{0}/\partial z^{2}\in{\mathscr{V}}_{n-2}$ and hence
orthogonal to $\mathbf{u}_{0}\in{\mathscr{V}}_{n}$. This explains why the
dissipation doesn’t vanish for a cylinder (inertial waves are not polynomials)
but does also for an ellipsoidal (inertial waves are polynomials in $x$, $y$
and $z$) - see §5 of Vantieghem (2014).

Figures 7 and 8 show the corrective viscous velocity defined in equation (22)
for the inertial mode (2,4,1). This is consistent with figure 6 in showing
that the diffusion-free boundary conditions lead to a smaller viscous
corrective flow than stress-free boundary conditions. In particular, the
stress-free boundary conditions give a corrective viscous velocity
$O(E^{1/2})$ smaller than the nonslip case everywhere due to the interior
viscous shear layers, while the diffusion-free solution is fully $O(E)$
smaller.

##  6 Inertial wave problem for a spherical shell

A spherical shell is an interesting singular geometry for the inertial wave
problem as smooth inertial modes do not exist in the inviscid limit except for
purely toroidal modes where the radial flow is zero everywhere. By adding a
small viscosity, inertial waves exist in the form of conical shear layers
spawned from the critical latitudes, though large scale structure may be
hidden beneath localized shear layers at certain frequencies (Lin & Ogilvie,
2021). Figure 9 shows the decay rate of an inertial wave with $m=2$ in a
spherical shell of inner-to-outer radius ratio of 0.5. Now there is no obvious
reduction in the corrective viscous velocity between the stress-free and
diffusion-free boundary conditions. This is because there is no benefit of
applying the latter if there is no solution of the inviscid problem. We would
expect a benefit for the toroidal modes but have not checked this.

Figure 9:  Decay rate of an inertial mode with $m=2$ and $\lambda=-1.4406$ in
a spherical shell of radius ratio 0.5.

##  7 Discussion

In this short note, we have discussed an apparently little-known type of
boundary condition that is a plausible alternative to the usual ‘stress-free’
or Neumann conditions used as the boundary conditions which do least harm to
the interior when diffusive processes are very weak. Originally used in
viscoelastic flow modelling, the objective here has been to show how these
‘diffusion-free” boundary conditions work for a couple of illustrative (ODE)
model problems and then demonstrate that they work, as expected, for a more
complicated PDE problem. The boundary layer correction is significantly weaker
for the diffusion-free boundary conditions except when the diffusion-free
problem doesn’t have a well-defined solution itself (the spherical shell).
Even in this latter scenario, diffusion-free does no worse than Neumann.

Hopefully, these examples highlight how useful these boundary conditions could
be in GAFD modelling. The boundary conditions can either be phrased as
‘turning off the diffusion at the wall’ or, equivalently for smooth solutions,
imposing that the ‘diffusion term vanishes’ which in its simplest form means
setting the second derivative to zero. The latter perspective is admittedly
unusual: no textbooks discuss setting the second derivative to vanish at a
boundary for a second order equation. Certainly for linear equations, this is
because, of course, this is really just asking the linear combination of lower
derivatives which define the second order derivative to vanish which is a
Robin condition. The point of the first model problem is, however, to show
that this can be generalised to the nonlinear setting where the exact
relationship is defined by turning diffusion off.

## Acknowledgements

RRK would like to thank Andrew Jackson and Greg Chini for indirectly
stimulating this note which was written for Andrew Soward’s 80th birthday
meeting held at Newcastle University in January 2024. The authors also thank
Toby Wood for bringing the existence of the Ventcel/Ventsell/Wentzell boundary
condition to our attention and Stephen K. Wilson for suggesting the ‘do least
harm’ moniker.

## Disclosure statement

No potential conflict of interest was reported by the authors.

## Funding

YL was supported by the National Natural Science Foundation of China (Nos.
12250012, 42142034).

## References

  * Bryan (1889) Bryan, G.H., The Waves on a Rotating Liquid Spheroid of Finite Ellipticity. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 1889, 180, 187–219.
  * Buffe (2017) Buffe, R. “Stabilization of the wave equation with Ventcel boundary condition” Journal de Mathematiques Pures et Appliquees, 108, 207-259, (2017)
  * Denk et al. (2021) Denk, R. and Kunze, M. and Ploss, D. “The bi-Laplacian with Wentzell boundary conditions on Lipschitz domains” Integral Equations and Operator Theory, 93, 13 (26 pages), (2021)
  * Feller (1952) Feller, W. “The parabolic differential equations and the asscoiated semi-groups of transformations” Ann. of Math. 55, 468-519, (1952)
  * Greenspan (1968) Greenspan, H. P. “The Theory of Rotating Fluids” Cambridge University Press, 1968.
  * He et al. (2022) He, J. and Favier, B. and Rieutrod, M. and Le Dizes, S. “Internal shear layers in librating spherical shells: the case of periodic characteristic paths” J. Fluid Mech. 939 A3 (30 pages), 2022.
  * He et al. (2023) He, J. and Favier, B. and Rieutord, M. and Le Dizes, S. “Internal shear layers in librating spherical shells: the case of attractors” J. Fluid Mech. submitted, 2023.
  * Hinz et al. (2018) Hinz, M. and Lancia, M. R. and Teplyaev, A. and Vernole, P. “Fractal snowflake domain diffusion with boundary and interior drifts” J. Mathematical Analysis and Applications, 2018, 672-693, 2018.
  * Le Dizes (2023) Le Dizes, S. “Critical slope singularities in rotating and stratified fluids”, preprint 2023.
  * Kelvin (1880) Kelvin, Lord. “Vibrations of a columnar vortex” Phil. Mag., 10, 155-168, 1880.
  * Kerswell (1993) Kerswell, R. R. “The Instability of precessing flow’ Geophys. Astrophys. Fluid Dynam. 72 107-144, 1993.
  * Kerswell & Barenghi (1995) Kerswell, R. R. and Barenghi, C. F. “Decay of inertial waves in a rotating circular cylinder’ J. Fluid Mech. 285 203-214, 1995.
  * Kerswell (1999) Kerswell, R. R. “Secondary instabilities in rapidly rotating fluids: inertial wave breakdown’ J. Fluid Mech. 382 283-306, 1999.
  * Lin & Ogilvie (2021) Lin, Y. and Ogilvie, G. I. “Resonant tidal responses in rotating fluid bodies: global modes hidden beneath localised wave beams’ Ap. J. Lett. 918 L21 (7 pages), 2021.
  * McEwan (1970) McEwan, A. D. Inertial oscillations in a rotating fluid cylinder” J. Fluid Mech. 40, 603-640, 1970.
  * Poincare (1910) Poincaré, H. “Sur la précession des corps déformables” Bull. Astronomique, 27, 321-356, 1910.
  * Reese et al. (2006) Reese, D., Lignières, F. and Rieutord, M., Acoustic oscillations of rapidly rotating polytropic stars. Astronomy & Astrophysics, 2006, 455, 621–637.
  * Sureshkumar & Beris (1995) Sureshkumar, R and Beris, Antony N “Effect of artificial stress diffusivity on the stability of numerical calculations ans the dynamics of time-dependent viscoelastic flows” J. Non-Newtonian Fluid Mech., 60, 53-80, 1995.
  * Sureshkumar et al. (1997) Sureshkumar, R. and Beris, A. N. and Handler, R. A., “Direct numerical simulation of the turbulent channel flow of a polymer solution” Phys. of Fluids, 9, 743-755, 1997.
  * Vantieghem (2014) Vantieghem, S. “Inertial modes in a rotating triaxial ellipsoid’ Proc. Roy. Soc. A 470 20140093, 2014.
  * Ventcel (1959) Ventcel, A. D. “On boundary conditions for multidimensional diffusion processes” Theor. Probab. Appl. IV, no 2, 164-177, (1959)
  * Vidal & Cebron (2023) Vidal, J. and Cebron, D. “Precession-driven flows in stress-free ellipsoids’ J. Fluid Mech. 954 A5, 2023.
  * Zhang & Liao (2017) Zhang, K. and Liao, X. “Theory and Modelling of Rotating Fluids” Cambridge University Press, 2017.
  * Zhang et al. (2001) Zhang, K. and Earnshaw, P. and Liao, X. and Busse, F. H. “On inertial waves in a rotating fluid sphere” J. Fluid Mech. 437, 103-119, 2001.
  * Zhang et al. (2004) Zhang, K. and Liao, X. and Earnshaw, P. “On inertial waves and oscillations in a rapidly rotating fluid spheroid” J. Fluid Mech. 504, 1-40, 2004.


# LUXE: A new experiment to study non-perturbative QED in $e^{-}$-laser and
$\gamma$-laser collisions.

Ruth Jacobs  on behalf of the LUXE collaboration  Deutsches Elektronen-
Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

###### Abstract

The LUXE experiment (Laser Und XFEL Experiment) is a new experiment in
planning at DESY Hamburg using the electron beam of the European XFEL. At
LUXE, the aim is to study collisions between a high-intensity optical laser
and up to $16.5\,$GeV electrons from the Eu.XFEL electron beam, or,
alternatively, high-energy secondary photons. The physics objectives of LUXE
are to measure processes of Quantum Electrodynamics (QED) at the strong-field
frontier, where QED is non-perturbative. This manifests itself in the creation
of physical electron-positron pairs from the QED vacuum. LUXE intends to
measure the positron production rate in a new physics regime at an
unprecedented laser intensity. Additionally, the high-intensity Compton photon
beam of LUXE can be used to search for physics beyond the Standard Model.

##  1 Strong-field QED

The theory of Quantum Electrodynamics (QED) has been tested with unrivalled
precision. Most experimental tests of QED to this day, however, consider QED
only in the perturbative regime. LUXE [1] aims to probe QED in a new regime of
strong fields, created in interactions between the $16.5\,$GeV Eu.XFEL
electron beam and a $40\,\textrm{TW}$ optical laser, as well as a secondary
high-energy gamma photon beam and the laser. The Lorentz boost of the
electrons in conjunction with the high-intensity laser create an
electromagnetic field strength above the Schwinger limit [2] (in case of an
electric field $E_{\textrm{cr}}=m^{2}_{e}c^{3}/(e\hbar)\approx 1.32\times
10^{18}\,\textrm{V/m}$) 111Here, $m_{e}$ denotes the electron mass, $c$ the
speed of light in vacuum, $e$ the electron charge and $\hbar$ the reduced
Planck constant., where QED becomes non-perturbative. Above the Schwinger
limit, the work by the field over one Compton wavelength ${\mkern
0.75mu\mathchar 22\relax\mkern-9.75mu\lambda}=\hbar/(m_{e}c)$ of the probe
electron is larger than the rest mass $m_{e}$. As a result, the QED vacuum
becomes polarized and processes to all orders of the electron-laser coupling
contribute to the scattering amplitude [3].

Two parameters govern the strong-field QED physics of LUXE: The classical
nonlinearity parameter $\xi=m_{e}c^{2}E/(\hbar\omega_{L}E_{\textrm{cr}})$,
where $\omega_{L}$ is the laser wavelength, $E$ is the electromagnetic field
strength, and $E_{\textrm{cr}}$ is the Schwinger critical field, quantifies
the coupling between the laser background field and the probe particle. The
probability of an n-th order electron-laser interaction process is therefore
proportional to $\xi^{2n}$. The quantum non-linearity parameter
$\chi=eE{\mkern 0.75mu\mathchar 22\relax\mkern-9.75mu\lambda}/(mc^{2})$
characterizes the field strength experienced by the probe particle in its rest
frame in relation to the Schwinger critical field, as well as the recoil
experienced by the probe particle emitting a photon. The two parameters are
related via the equation $\chi=\xi\eta$, where
$\eta=\gamma\hbar\omega_{L}(1+\cos\theta)/(m_{e}c^{2})$ is the energy
parameter, $\gamma$ is the relativistic Lorentz factor, and $\theta$ is the
electron-laser crossing angle ($\theta\approx 17.2^{\circ}$ in LUXE).

The main processes of strong-field QED probed by LUXE are non-linear Compton
scattering (see fig. 1 (left)) and non-linear Breit-Wheeler pair creation (see
fig. 1 (right)). In non-linear Compton scattering, the probe electron absorbs
multiple laser photons and emits a single high-energy photon. The distinct
feature of non-linear Compton scattering is the displacement of the Compton
edge as function of the laser intensity parameter $\xi$, due to the fact that
the probe electron acquires a larger effective mass
$m_{e}^{*}=m_{e}\sqrt{1+\xi^{2}}$ in the laser field. Figure 2 (left) shows
the Compton electron energy spectrum for different values of $\xi$ simulated
with the Ptarmigan strong-field QED generator package [5].

Breit-Wheeler pair production is the field-induced creation of a physical
electron-positron pair from the QED vacuum. In this process, a high-energy
photon, produced in the LUXE $e^{-}$-laser collision mode via non-linear
Compton scattering, or, in the LUXE $\gamma$-laser mode, stemming from a
secondary high-energy photon beam, absorbs multiple laser photons and produces
an electron-positron pair. This process has no classical equivalent, unlike
the Compton process, therefore it is a direct probe of non-perturbative QED.
The positron production rate, $\Gamma_{\textrm{BW}}$ is one of the main
quantities of interest for LUXE. Figure 2 (right) shows the behaviour of
$\Gamma_{\textrm{BW}}$ as a function of the laser intensity $\xi$. For $\xi\ll
1$ the positron rate $\Gamma_{\textrm{BW}}$ follows a power law
$\Gamma_{\textrm{BW}}\sim\xi^{2n}$. In the regime $\xi\gg 1$, a departure from
the power-law characteristic occurs, where, for $\chi\ll 1$,
$\Gamma_{\textrm{BW}}$ scales according to:

$\Gamma_{\textrm{BW}}\propto\exp\left[-\frac{3}{8}\frac{1}{\omega_{\gamma}(1+\cos\theta)}\frac{E_{\textrm{cr}}}{E_{\textrm{L}}}\right].$
(1)

In the fully non-perturbative regime, where $\xi\gg 1$ and $\chi\gg 1$,
radiative quantum corrections to all orders contribute to
$\Gamma_{\textrm{BW}}$. An equivalent to eq. 1 in the case where $\chi\gtrsim
1$ can be found in the following article [4].

Figure 1: Feynman diagrams of strong-field QED processes at LUXE, showing non-
linear Compton scattering (left) and non-linear Breit-Wheeler pair production
(right).

Figure 2: Left: Simulated Compton electron energy spectrum for one electron-
laser bunch crossing as a function of the laser intensity $\xi$ in the LUXE
$e^{-}$-laser mode. Right: Breit-Wheeler positron rate as a function of $\xi$
in the LUXE $\gamma$-laser mode. Uncertainty corresponds to one week of data-
taking.

##  2 LUXE experimental setup

Figure 3 shows the LUXE experimental setup for the $e^{-}$-laser colision mode
(left) and the $\gamma$-laser collision mode (right). In the $e^{-}$-laser
mode, the $16.5\,\textrm{GeV}$ electron beam directly travels to the
interaction point where it collides with the high-intensity laser. In the
$\gamma$-laser mode, a secondary high-energy photon beam is created by
impinging the electron beam on a Tungsten Bremsstrahlung target, or,
alternatively,by colliding the beam electrons with a low-intensity laser pulse
to produce Compton photons with a narrow energy bandwidth.

The laser system foreseen for LUXE is a commercially sourced titanium-sapphire
($\lambda=800\,\textrm{nm}$) $40\,\textrm{TW}$ pulsed laser system (LUXE
phase-0), which is upgradeable to $350\,\textrm{TW}$ (LUXE phase-1). The
parameter space in terms of the nonlinearity parameters $\xi$ and $\chi$ that
is accessible with the phase-0 (phase-1) laser setup is $\xi<7.9$ ($\xi<23.6$)
and $\chi<1.5$ ($\chi<4.45$). One of the main challenges of LUXE is the
required shot-to-shot stability of the laser ($<1\%$ variation, $<5\%$
intensity uncertainty), which is ensured by a dedicated suite of laser
diagnostics.

The particle detectors in LUXE are used to detect electrons, positrons and
photons and to measure their energy based on the transverse displacement in a
dipole spectrometer. Since the particle rates per bunch-crossing (BX) in LUXE
are vastly different, depending on the run mode, location in the experimental
setup and laser intensity, dedicated technologies are used for each system.
For the positron detection ($10^{-3}<N_{e^{+}}/BX<10^{4}$) a 4-layer silicon
tracker in combination with a high-granularity electromagnetic calorimeter is
used. For the Compton electron detection system the challenge is to
reconstruct extremely high electron rates ($10^{3}<N_{e^{-}}/BX<10^{8}$). A
finely segmented air-filled Cherenkov detector in combination with a
scintillator screen read out by an optical camera system provides a robust
solution. Finally, the photon energy spectrum is measured by a combination of
three complementary detectors in the forward region of the experiment.
Firstly, the gamma spectrometer partly reconverts the photon beam on a target
into electron-positron pairs, which are subsequently analysed in a dipole
spectrometer using a scintillator and camera system similar to the one used
for the Compton-scattered electrons. Secondly, the gamma beam profiler studies
the beam distribution in the transverse plane using sapphire strip detectors,
in order to precisely determine the laser intensity $\xi$. Finally, a lead-
glass calorimeter monitors backscattered particles from impinging the photons
on the final dump, thereby measuring the total photon flux. In general, at
least two complementary detector technologies are foreseen for each location
in LUXE, to enable cross-calibration and the reduction of systematic
uncertainties.

In addition to the LUXE strong-field QED program, the experimental setup can
be extended to probe physics beyond the Standard Model (BSM) coupling to
photons [6]. Rare BSM particles, such as axion-like particles (ALPs) could be
produced in the interaction of the LUXE high-intensity Compton photon beam
with the photon beam dump material, for example via Primakov production. A
photon detector placed in a well-shielded area behind the beam dump can be
used to search for ALPs decaying to two photons in the volume between the dump
and the detector. The expected sensitivity of the LUXE extended setup is
comparable to ongoing and planned experiments searching for ALPs.

Figure 3: Left: Schematic of the LUXE experimental setup in $e^{-}$-laser
mode. Right: LUXE experimental setup in $\gamma$-laser mode.

##  3 Conclusion

The aim of LUXE is to probe QED in a new regime of strong-fields, by studying
collisions between the EU.XFEL electron beam, or a high-energy secondary gamma
photon beam, with a high-intensity optical laser. Running LUXE as a collision
experiment in continuous data-taking mode will enable precision measurements
of strong-field QED processes, such as non-linear Compton scattering and
Breit-Wheeler pair production. The laser system and particle detectors in LUXE
are custom-designed in order to meet the physics goals. LUXE will likely be
the first experiment to take precision measurements in a regime of QED never
before explored in clean laboratory conditions and to study high-intensity
laser collisions with real high-energy gamma photons.

## Acknowledgments

This work was in part funded by the Deutsche Forschungsgemeinschaft under
Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306.

This work has benefited from computing services provided by the German
National Analysis Facility (NAF).

## References

## References

  * [1] LUXE collaboration, EPJ S. T. 230, 11 (2021).
  * [2] J. Schwinger, Phys. Rev. 82, 664 (1951).
  * [3] A. Fedotov et al., RIKEN-iTHEMS-Report-22, https://arxiv.org/abs/2203.00019, 2022.
  * [4] T. G. Blackburn, B. King, EPJ C 82, 44 (2022)
  * [5] T.G. Blackburn, https://github.com/tgblackburn/ptarmigan.
  * [6] Bai, Z. et al., https://arxiv.org/abs/2107.13554, 2021.


Note that in the breadth-first-search and blocking-flow algorithms in <ref>, there is an $O(k)$ overhead where we have to spend $O(k)$ time iterating through the $k$ binary search trees in order to explore the out-neighbors of the $O(kr)$ elements.
Our goal in this section is thus to show that it is possible to further “sparsify” the exchange graphs to contain essentially only a basis, hence reducing its size from $\Theta(kr)$ to $O(r)$.
We start with a slight modification to the BFS <ref> which reduces the running time by a factor of $O(k)$. The idea is that if we visit an element $u$ in the BFS which does not increase the rank of all visited elements so far, we can skip searching out-edges from $u$. Indeed, if $(u,v)$ is an edge of the exchange graph, then there must have been some element $u'$ visited earlier in the BFS which also has the edge $(u',v)$.

BFS in a $k$-fold union exchange graph

$S \subseteq U$ with partition $S_1, \ldots, S_k$ of independent sets and a basis $B$ of $U \setminus S$
The $(s, v)$-distance $d(v)$ in $H(S)$ for each $v \in S \cup \{t\}$

$\mathsf{queue} \gets B$ and $R \gets \emptyset$
$d(v) \gets \infty$ for each $v \in S \cup \{t\}$, and $d(v) \gets 1$ for each $v \in B$
$\cT_i \gets \textsc{Initialize}(\cM, S_i, S_i)$ (<ref> with $\beta = 1$)
$\mathsf{queue} \neq \emptyset$
$u \gets \mathsf{queue}.\textsc{Pop}()$
$R + u \in \cI$
$i \in \{1, 2, \ldots, k\}$
$v := \cT_i.\textsc{Find}(u) \neq \bot$
$d(v) \gets d(u) + 1$ and $\mathsf{queue}.\textsc{Push}(v)$

$S_i + u \in \cI$ and $d(t) = \infty$
$d(t) \gets d(u) + 1$

$R \gets R + u$

return $d(v)$ for each $v \in S \cup \{t\}$.

Given $S \in \cI_\text{part} \cap \hat{\cI}$ and a basis $B$ of $U \setminus S$, it takes $\tO(kr)$ time to construct the distance layers $L_2, \ldots, L_{d_t - 1}$ of $H(S)$.

The algorithm is presented as <ref>.
It performs a BFS from a basis $B$ of the first layer and only explores out-edges from the first basis $R$ it found.
It takes
(i) $\tO(|S|)$ time to construct the $\cT_i$'s,
(ii) $\tO(1)$ time to discover each of the $O(kr)$ element, and
(iii) an additional $\tO(k \cdot |R|)$ time to iterate through all $k$ binary search trees $\cT_i$'s for each $u \in R$.
The total running time is thus bounded by $\tO(kr)$.

We have shown in <ref> that it is feasible to replace $U \setminus S$ with simply $B$.
It remains to show that exploring only the out-neighbors of $u \in R$ does not affect the correctness.
Consider a $v \in S \setminus R$ (we know that $B \subseteq R$ so it suffices to consider elements in $S$) with an out-neighbor $x \in S_i$, i.e., $S_i - x + v \in \cI$.
It then follows that $\rank((S_i - x) + R_v) = \rank((S_i - x) + (R_v + v)) \geq \rank(S_i)$ by <ref> and <ref>, where $R_v$ is the set $R$ when $v$ is popped out of the queue (in other words, $R_v + v \not\in \cI$).
This implies that there is a $u \in R_v$ which is visited before $v$ that also has out-neighbor $x$.
The modification is therefore correct.

Our blocking-flow algorithm for $k$-fold matroid union is presented as <ref>.
It's essentially a specialization of <ref> to the case where all the $k$ matroids are the same, except that we skip exploring the out-neighbors of $a_{\ell}$ and remove it directly if it is “spanned” by the previous layers and the set $R_{\ell} \subseteq L_{\ell}$ of elements that are not on any augmenting path of length $d_t$.
With this optimization, we obtain the following lemma analogous to <ref>.

Blocking flow in a $k$-fold union exchange graph

$S \subseteq U$ which partitions into $S_1, \ldots, S_k$ of independent sets and a dynamic-basis data structure $\cD$ of $U \setminus S$
$S^\prime \in \cI_{\text{part}} \cap \cI_k$ with $d_{H(S^\prime)}(s, t) > d_{H(S)}(s, t)$
$\cD$ maintains a basis of $U \setminus S^\prime$ at the end of the algorithm

Build the distance layers $L_2, \ldots, L_{d_t - 1}$ of $H(S)$ with <ref>
$L_0 \gets \{s\}$ and $L_{d_t} \gets \{t\}$
$B \gets$ the basis maintained by $\cD$ and $L_1 \gets B$
$A_\ell \gets L_\ell$ for each $0 \leq \ell \leq d_t$
$\cT_{\ell}^{(i)} \gets \textsc{Initialize}(\cM_i, S_i, Q_{S_i}, A_{\ell} \cap S_i)$ for each $2 \leq \ell < d_t$ and $1 \leq i \leq k$ (<ref> with $\beta = \sqrt{r} / d_t$)
$D_{\ell} \gets \emptyset$ for each $1 \leq \ell < d_t$
$R_{\ell} \gets \emptyset$ for each $2 \leq \ell < d_t$
$\ell \gets 0$ and $a_0 \gets s$
$\ell \geq 0$
$\ell < d_t$
$A_{\ell} = \emptyset$

$\ell \geq 2$ and $\rank(L_1 \cup \cdots \cup L_{\ell - 1} \cup R_{\ell} \cup \{a_{\ell}\}) = \rank(L_1 \cup \cdots \cup L_{\ell - 1} \cup R_{\ell})$
$A_{\ell} \gets A_{\ell} - a_{\ell}$ and continue

$\ell > 0$
Find an $a_{\ell + 1} := \cT_{\ell + 1}^{(i)}.\textsc{Find}(a_\ell) \neq \bot$ for some $1 \leq i \leq k$

$a_{\ell + 1} \gets$ an arbitrary element in $A_1$

such an $a_{\ell + 1}$ does not exist
$\ell \geq 2$
$R_{\ell} \gets R_{\ell} + a_{\ell}$ and $\cT_{\ell}^{(j)}.\textsc{Delete}(a_{\ell})$ where $a_{\ell} \in S_j$

$A_\ell \gets A_\ell - a_\ell$ and $\ell \gets \ell - 1$

$\ell \gets \ell + 1$

Found augmenting path $a_1, a_2, \ldots a_\ell$
$B \gets B - a_1$, $A_1 \gets A_1 - a_1$, and $D_1 \gets D_1 + a_1$
$\cD.\textsc{Delete}(a_1)$ returns a replacement $x$
$B_i \gets B_i + x$  and $A_i \gets A_i + x$

$i \in \{2, \ldots, d_t - 1\}$
$D_{i} \gets D_{i} + a_{i}$ and $A_{i} \gets A_{i} - a_{i}$
$\cT_{i}^{(j)}.\textsc{Delete}(a_i)$ and $\cT_{i}^{(j)}.\textsc{Update}(\{a_{i - 1}, a_{i}\})$ where $a_i \in S_j$

Augment $S$ along $P = (s, a_1, \ldots, a_{d_t - 1}, t)$
$\ell \gets 0$

return $S$

Given an $S \in \cI_{\text{part}} \cap \hat{\cI}$ with $d_{H(S)}(s, t) = d_t$ together with a data structure $\cD$ of <ref> that maintains a basis of $U \setminus S$, it takes

\begin{equation}
\tO\left(\underbrace{kr}_{\ref{item:term1}} + \underbrace{\left(|S^\prime| - |S|\right) \cdot d_t\sqrt{r}}_{\ref{item:term2}} + \underbrace{\left((|S^\prime| - |S|) \cdot d_t + r\right) \cdot k}_{\ref{item:term3}} + \underbrace{\left(kr + (|S^\prime| - |S|)\right) \cdot \frac{\sqrt{r}}{d_t}}_{\ref{item:term4}}\right)
\label{eq:complexity}
\end{equation}

time to obtain an $S^\prime \in \cI_{\text{part}} \cap \hat{\cI}$ with $d_{H}(S^\prime)(s, t) > d_t$, with an additional guarantee that $\cD$ now maintains a basis of $U \setminus S^\prime$.

We need the following observation to bound the running time of <ref>.

In <ref>, it holds that $B \cup R_2 \cup R_{3} \cup \cdots \cup R_{d_t - 1} \in \cI$.

We analyze the running time of <ref> first.
In particular, there are four terms in <ref> which come from the following.

*  $\tO(kr)$: It takes $\tO(kr)$ time to compute the distance layers using <ref> and initialize all the binary search trees $\cT_{\ell}^{(i)}$'s. Computing the rank of $L_1 \cup \cdots \cup L_{\ell - 1} \cup R_{\ell}$ also takes $\tO(kr)$ time in total since we can pre-compute query-sets of the form $L_1 \cup \cdots \cup L_{k}$ for each $k$ in $\tO(kr)$ time, and each insertion to $R_{\ell}$ takes $\tO(1)$ time.

*  $\tO\left(\left(|S^\prime| - |S|\right) \cdot d_t\sqrt{r}\right)$: For each of the $O(|S^\prime| - |S|)$ augmentations, it takes $\tO(r \cdot \frac{d_t}{\sqrt{r}})$ time to update the binary search trees.

*  $\tO(\left((|S^\prime| - |S|) \cdot d_t + r\right) \cdot k)$: The number of elements whose out-edges are explored is bounded by $O\left((|S^\prime| - |S|) \cdot d_t + r\right)$. This is because for each such element $u$, either $u$ is included in an augmenting path of length $d_t$, or $u$ is removed in Line <ref>.
There are $O((|S^\prime| - |S|) \cdot d_t)$ such $u$'s in the augmenting paths.
For $u$ removed in Line <ref>, if $\ell = 1$, then the number of such $u$'s is $O(|S^\prime| - |S| + r)$ because there are initially $O(r)$ elements in $A_1$, and we add at most one to it every augmentation.
If $\ell \geq 2$, then we insert it into $R_{\ell}$, and by Line <ref>, the rank of $L_1 \cup \cdots \cup L_{\ell - 1} \cup R_{\ell}$ increases after including $u$ into $R_{\ell}$.
By <ref>, the number of such $u$'s is bounded by $O(r)$.
The term then comes from spending $O(k)$ time iterating through the $k$ binary search trees for each of the $O\left((|S^\prime| - |S|) \cdot d_t + r\right)$ elements whose out-neighbors are explored.

*  $\tO(\left(kr + (|S^\prime| - |S|)\right) \cdot \frac{\sqrt{r}}{d_t})$:  The number of elements that are once in some $A_{\ell}$ is bounded by $O(kr + |S^\prime| - |S|)$. Initially, there are $O(kr)$ elements ($A_1$ plus all the $A_\ell$ for $\ell \geq 2$), and each augmentation adds at most one element to $A_1$.
Each of these elements is discovered by $\cT_{\ell}^{(i)}.\textsc{Find}(\cdot)$ at most once, and thus we can charge the $\tO(\frac{\sqrt{r}}{d_t})$ cost to it, resulting in the fourth term of <ref>.

Note that for each element whose out-neighbors are explored, any failed attempt of $\cT_{\ell}^{(i)}.\textsc{Find}(\cdot)$ costs only $\tO(1)$ instead of $\tO(\frac{\sqrt{r}}{d_t})$ according to <ref>.
The $\tO(\frac{\sqrt{r}}{d_t})$ cost of a successful search is charged to term <ref> instead of <ref>.

As for correctness, it suffices to show that each of the $a_{\ell}$ removed from $A_{\ell}$ because it is spanned by $L_1 \cup \cdots \cup L_{\ell - 1} \cup R_{\ell}$ in Line <ref> is not in any augmenting path of length $d_t$.
Consider its out-neighbor $a_{\ell + 1}$ with respect to the current $S$, and we would like to argue that $a_{\ell + 1}$ is not on any augmenting path of length $d_t$ anymore.
This is because we have already explored all the out-neighbors of elements in $R_{\ell}$.
Since $a_{\ell} \in \sspan(L_1 \cup \cdots \cup L_{\ell - 1} \cup R_{\ell})$, by <ref>, there must exist some $u \in L_1 \cup \cdots \cup L_{\ell - 1} \cup R_{\ell}$ with a directed edge $(u, a_{\ell + 1})$.
We consider two cases:

* $u\in R_{\ell}$. This means that we have already explored $a_{\ell+1}$, as we finished exploring all out-neighbors of $u$ already.

* $u\in L_{1}\cup \cdots \cup L_{\ell-1}$. We know that by <ref>, both $d_{H(S)}(s,v)$ and
$d_{H(S)}(v,t)$ can only increase after augmentations for all elements $v$.
Hence $a_{\ell+1}$ cannot be part of an augmenting path of length $d_t$ anymore, since if it was its distance to $t$ must be $d-(\ell+1)$, but then the distance from $u$ to $t$ must be at most $d-\ell$ (which is smaller than its initial distance to $t$ at the beginning of the phase).

As a result, all of $u$'s out-neighbors have either already been explored or do not belong to any augmenting path of length $d_t$.
This implies that $u$ is not on any such path either, and thus it's correct to skip and remove it from $A_{\ell}$.
This concludes the proof of <ref>.

<ref> now follows from analyzing the total running time of $O(\sqrt{\min(n, kr)})$ runs of <ref>.

We initialize the dynamic-basis data structure $\cD$ of <ref> on $U$ in $\tO(n)$ time.
Let $p = \min(n, kr)$ be the rank of the $k$-fold matroid union.
Using $\cD$, we then run $O(\sqrt{p})$ iterations of <ref> until $d_{H(S)}(s, t) > \sqrt{p}$.
Summing the first two terms of <ref> over these $O(\sqrt{p})$ iterations gives (recall that <ref> guarantees that $\sum_{d = 1}^{\sqrt{p}}d \cdot (|S_d| - |S_{d - 1}|) = \tO(p)$)

\[
\tO\left(kr\sqrt{p} + \sqrt{r} \cdot \sum_{d = 1}^{\sqrt{p}}d \cdot \left(|S_d| - |S_{d - 1}|\right)\right) = \tO\left(kr\sqrt{p}\right)
\]

since $p\sqrt{r} \leq kr\sqrt{p}$.
The third term of <ref> contributes a total running time of

\[
\tO\left(\left(\sum_{d = 1}^{\sqrt{p}}dk \cdot \left(|S_d| - |S_{d - 1}|\right)\right) + kr\sqrt{p}\right) = \tO\left(kr\sqrt{p} + kp\right),
\]

while the fourth term of <ref> sums up to

\[
\tO\left(\left(\sum_{d = 1}^{\sqrt{p}}kr\frac{\sqrt{r}}{d}\right) + kr\sqrt{r}\right) = \tO\left(kr\sqrt{r}\right).
\]

We finish the algorithm by finding the remaining $O(\sqrt{p})$ augmenting paths one at a time with <ref> in a total of $\tO(kr\sqrt{p})$ time.
The $k$-fold matroid union algorithm thus indeed runs in $\tO\left(n + kr\sqrt{\min(n, kr)} + k\min(n, kr)\right)$ time, concluding the proof of <ref>.

§ DYNAMIC ORACLES FOR SPECIFIC MATROIDS & APPLICATIONS

In this appendix, we show how to leverage known dynamic algorithms to implement the dynamic rank oracle (<ref>) efficiently for many important matroids.
What we need are data structures that can maintain the rank of a set dynamically under insertions and deletions in worst-case update time (converting a worst-case data structure to fully-persistent can be done by the standard technique of [29, 27], paying an overhead of $O(\log{n})$). Additionally, note that the data structures do not need to work against an adaptive adversary since we only ever use the rank of the queried sets, which is not affected by internal randomness.

In particular, for partition, graphic, bicircular, convex transversal, and simple job scheduling matroids it is possible to maintain the rank with polylog$(n)$ update-time, and for linear matroids in $O(n^{1.529})$ update-time.

Together with our matroid intersection (<ref>) and matroid union (<ref>) algorithms, this leads to a black-box approach to solving many different problems. In fact, we can solve matroid intersection and union on any combination of the above matroids, leading to improved or matching running times for many problems (see the introduction <ref> with <ref> for a more thorough discussion). For completeness, we define these problems in <ref>.
The same algorithms are powerful enough to also solve new problems which have not been studied before.

Example Application: Tree Scheduling (or Maximum Forest with Deadlines). We give an example of a reasonably natural combinatorial scheduling problem, which—to our knowledge—has not been studied before. Suppose we are given a graph $G = (V,E)$ where each edge $e\in E$ has two numbers associated with it: a release date $\ell_e$ and a deadline $r_e$. Consider the problem where we want to for each day pick exactly one edge (say, to build/construct), but we have constraints that edge $e$ can only be picked between days $\ell_e$ and $r_e$. Now the task is to come up with a scheduling plan to build a spanning tree of the graph, if possible.

This problem is exactly a matroid intersection problem between a graphic matroid and a convex transversal matroid. Hence, by a black-box reduction, we know that we can solve this problem in $\tO(|E|\sqrt{|V|})$ time.

§.§ Partition Matroids

In a partition matroid $\cM = (\ground,\cI)$, each element $u\in \ground$ is assigned a color $c_u$. We are also, for each color $c$, given a non-negative integer $d_c$, and we define a set of elements $S\subseteq \ground$ to be independent if for each color $c$, $S$ includes at most $d_c$ elements of this color.
Implementing the dynamic oracle for the partition matroid is easy:

One can maintain the rank of a partition matroid in $O(1)$-update time.

For each color $c$ we maintain a counter $x_c$ of how many elements we have of color $c$. We also maintain $r = \sum \min_c(x_c, d_c)$, which is the rank of the current set.

Bipartite matching can be modeled as a matroid intersection problem of two partition matroids. So our matroid intersection algorithm together with the above lemma match (up to poly-logarithmic factors induced by fully-persistence)—in a black-box fashion—the $O(|E|\sqrt{|V|})$-time bound of the best combinatorial algorithm for bipartite matching [49].

§.§ Graphic and Bicircular Matroids

Given a graph $G = (V,E)$, the graphic and bicircular matroids are matroids capturing the connectivity structure of the graph.

Graphic Matroid.
In the graphic matroids $\cM = (E,\cI)$ a subset of edges $E'\subseteq E$ are independent if and only if they do not contain a cycle.
We use the following result to implement the dynamic oracle for this matroid.

There is a data structure that maintains an initially-empty graph $G = (V, E)$ and supports insertion/deletion of edges $e$ into/from $E$ in worst case $O(\log^4{|V|})$ time and query of the connectivity between $u$ and $v$ in worst case $O(\log{|V|}/\log{\log{|V|}})$ time.
The data structure works with high probability against an oblivious adversary.

With a simple and standard extension, we can maintain the number of connected components as well, and hence also the rank (since $\rank(E') = |V|- \#\text{connected components in }G[E']$).

There is a data structure that maintains an initially-empty graph $G = (V, E)$ and supports insertion/deletion of $e$ into/from $E$ in worst-case $O(\log^4{|V|})$ time.
After each operation, the data structure also returns the number of connected components in $G$.
The data structure works with high probability against an oblivious adversary.

We maintain the data structure $\cC$ of <ref> and a counter $c := |V|$ representing the number of connected components.
For insertion of $e = (u, v)$, we first query the connectivity of $u$ and $v$ before inserting $e$ into $\cC$.
If they are not connected before the insertion, decrease $c$ by one.
For deletion of $e = (u, v)$, after deleting $e$ from $\cC$, we check if $u$ and $v$ are still connected.
If not, then we increase $c$ by one.

Bicircular Matroid.
In the bicircular matroid $\cM = (E,\cI)$, a subset of edges $E'\subseteq E$ are independent if and only if each connected component in $G[E']$ has at most one cycle. Similar to the graphic matroid, dynamic connectivity algorithms can be used to implement the dynamic rank oracle for bicircular matroids too.

There is a data structure that maintains an initially-empty graph $G = (V, E)$ and supports insertion/deletion of $e$ into/from $E$ in worst-case $O(\log^4{|V|})$ time.
After each operation, the data structure also returns the rank of $E$ in the bicircular matroid.
The data structure works with high probability against an oblivious adversary.

The dynamic connectivity data structure of
[52, 40] (<ref>) can be adapted to also keep track of the number of edges and vertices in each connected component. Using this, the data structure can, for each connected component $c$ keep track of a number $x_c$ as the minimum of the number of edges in this component and the number of vertices in this component. Then the rank of the bicircular matroid is just the sum of $x_c$ (as in an independent set each component is either a tree or a tree with an extra edge). In each update two components can merge, a component can be split up into two, or the edge-count of a component may simply change.

The above dynamic connectivity data structures are randomized.
There are also deterministic connectivity data structures, but with slightly less efficient sub-polynomial $|V|^{o(1)}$ update time [19].

§.§ Convex Transversal and Scheduling Matroids

Convex transversal and  scheduling matroids are special cases of the transversal matroid, with applications in scheduling algorithms.

A transversal matroid with respect to a bipartite graph $G = (L, R, E)$ is defined over the ground set $L$, where each $S \subseteq L$ is independent if and only if there is a perfect matching in $G$ between $S$ and a subset of $R$.

A bipartite graph $G = (L, R, E)$ is convex if $R$ has a linear order $R = \{r_1, r_2, \ldots, r_n\}$ and each $\ell \in L$ corresponds to an interval $1 \leq s(\ell) \leq t(\ell) \leq n$ such that $(\ell, r_i) \in E$ if and only if $s(\ell) \leq i \leq t(\ell)$, i.e., the neighbors of each $\ell$ form an interval.

A convex transversal matroid is a transversal matroid with respect to a convex bipartite graph.
A simple job scheduling matroid is a special case of convex transversal matroids in which $s(\ell) = 1$ for each $\ell \in L$.

One intuitive way to think about the simple job scheduling matroid is that there is a machine capable of finishing one job per day.
The ground set of the matroid consists of $n$ jobs, where the $i$-th jobs must be done before its deadline $d_i$.
A subset of jobs forms an independent set if it's possible to schedule these jobs on the machine so that every job is finished before its deadline.

There is a data structure which, given a convex bipartite graph $G = (L, R, E)$, after $\tO(|L|+|R|)$ initialization, maintains the size of the maximum matching of $G[A \cup R]$ where $A \subseteq L$ is a dynamically changing subset of $L$ that is initially empty.
The data structure supports insertion/deletion of an $x \in L$ to/from $A$ in worst-case $O(\log^{2}(|L|+|R|))$ update time.

The exact data structure presented in [13] is different from the stated one.
In particular, they support insertion/deletion of an unknown job, i.e., we do not know beforehand what the starting date and deadline of the job are, nor do we know its relative position among the current set of jobs.
As a result, they used a rebalancing-based or rebuilding-based binary search tree [68, 10, 11], resulting in their amortized bound.
For our use case, all the possible jobs are known and we are just activating/deactivating them, hence a static binary tree with a worst-case guarantee over these jobs suffices.

§.§ Linear Matroid

In a linear matroid $\cM = (\ground, \cI)$, $\ground$ is a set of $n$ vectors (of dimension $r$) in some vector space and the notion of independence is just that of linear independence. The dynamic algorithm to maintain the rank of a matrix of [81] can be used without modification as the dynamic oracle.

There is a data structure which, given an $n \times n$ matrix $M$, maintains the rank of $M$ under row updates in worst-case $O(n^{1.529})$ update time.

§.§ Problems

For completeness, here we define the problems we discuss in the introduction, and why they reduce to matroid union or intersection.

$k$-Forest. In this problem we are given a graph $G = (V,E)$ and asked to find $k$ edge-disjoint forests of the graph, of the maximum total size.
It can be modeled as the $k$-fold matroids union over the graphic matroid of $G$.

$k$-Disjoint Spanning Trees. This problem is a special case of the above $k$-forest problem where we ask to find $k$ edge-disjoint spanning trees of the graph. Clearly, if such exists, the $k$-forest problem will find them.

$k$-Pseudoforest. Similar to above, in this problem we are given a graph $G = (V,E)$ and asked to find $k$ edge-disjoint pseudoforests of the graph, of the maximum total size. A pseudoforest is an undirected graph in which every component has at most one cycle.
The problem can be modeled as the $k$-fold matroids union over the bicircular matroid of $G$.

$(f,p)$-Mixed Forest-Pseudoforest. Again, we are given a graph $G = (V,E)$ and asked to find $f$ forests and $p$ pseudoforest (all edge-disjoint), of the maximum total size.
The problem can be modeled as the matroids union over $f$ graphic matroids and $p$ bicircular matroids.

Tree Packing. In the tree packing problem, we are given a graph $G = (V,E)$ and are asked to find the maximum $k$ such that we can find $k$-disjoint spanning trees in the graph. This number $k$ is sometimes called the tree-pack-number or strength of the graph. The problem can be solved with the $k$-disjoint spanning trees problem, by binary searching for $k$ in the range $[0,|E|/(|V|-1)]$, and is an example of a matroid packing problem.

Arboricity and Pseudoarboricity. The arboricity (respectively pseudoarboricity) of a graph $G = (V,E)$ is the least integer $k$ such that we can partition the edges into $k$ edge-disjoint forests (respectively pseudoforests).
This can be solved with the $k$-forest (respectively $k$-pseudoforest) problem with a binary search over $k$. It is well known that for a simple graph the (pseudo-)arboricity is at most $\sqrt{|E|}$, so we need only search for $k$ in the range $[0,\sqrt{|E|}]$. The problems are examples of matroid covering problems.

Shannon Switching Game. The Shannon switching game is a game played on a graph $G = (V,E)$, between two players “Short” and “Cut”. They alternate turns with Short playing first, and all edges are initially colored white. On Short's turn, he may color an edge of the graph black. On Cut's turn, he picks a remaining non-black edge and removes it from the graph. Short wins if he connects the full graph with black edges, and Cut if he manages to disconnect the graph. It can be shown that Short wins if and only if there exists two disjoint spanning trees in the graph (and these two spanning trees describes a winning strategy for Short). Hence solving this game is a special case of the $k$-disjoint spanning tree problem with $k = 2$.

Graph $k$-Irreducibility.
A (multi-)graph $G = (V,E)$ is called $k$-irreducible ([83]) if and only if $|E| = k(|V|-1)$ and for any vertex-induced nonempty, proper subgraph $G[V']$ it holds that $|E(G[V'])| < k(|V'|-1)$.
The motivation behind this definition comes from the rigidity of bar-and-body frameworks. A bar-and-body framework where rigid bars are attached to rigid bodied with joints (represented by the graph $G$).
Then any stress put on a $k$-irreducible structure will propagate to all the bars (i.e. edges).
[45] show how one can decide if a graph is $k$-irreducible by first determining if its edges can be partitioned into $k$ edge-disjoint trees, and then performing an additional $\tO(k|V|)$ work.

Bipartite Matching. In the bipartite matching problem, we are given a bipartite graph $G = (L\cup R,E)$, and the goal is to find a matching (a set of edges which share no vertices) of maximum size. Bipartite matching can be modeled as a matroid intersection problem over two partition matroids $M_L = (E,\cI_L)$ and $M_R = (E,\cI_R)$. $M_L$ specifies that no two edges share the same vertex on the left $L$ (and $M_R$ is defined similarly on the right  set of vertices $R$).

Colorful Spanning Tree. In this problem[sometimes also called rainbow spanning tree.], we are given a graph $G = (V,E)$ together with colors on the edges $c:E \to \bbZ$. We are tasked to find a spanning tree of $G$ such that no two edges in our spanning tree have the same color. This problem can be modeled by the matroid intersection of the graphic matroid of $G$ (ensuring we pick a forest),
and a partition matroid of the coloring $c$ (ensuring that we pick no duplicate colors). We also note that this problem is more difficult than bipartite matching since any bipartite matching instance can be converted to a colorful spanning tree instance on a star-multi-graph.

Graphic Matroid Intersection. In graphic matroid intersection we are given two graphs $G_1 = (V_1,E_1)$ and $G_2 = (V_2,E_2)$ and a bijection of the edges $\phi : E_1 \to E_2$. The task is to find a forest in $G_1$ of the maximum size, which also maps to a forest in $G_2$. By definition, this is a matroid intersection problem over two graphic matroids. Again, this problem is a further generalization of the colorful spanning tree problem.

Convex Transversal and Simple Job Scheduling Matroid Intersection. In these problems, we are given a set of unit-size jobs $V$, where each job $v$ has two release times $\ell_1(v)$, $\ell_2(v) \ge 1$ (in simple job scheduling $\ell(v) = 1$) and two deadlines $r_1(v)$,$r_2(v)\le \mu$. The task is to find a set of jobs $S$ of the maximum size such that they can be scheduled on two machines as follows: each job needs to be scheduled at both machines, and at machine $i$ it must be scheduled at time $t \in [\ell_i(v), r_i(v)]$.

Linear Matroid Intersection. In this problem, we are given two $n\times r$ matrices $M_1$ and $M_2$ over some field.
The task is to find a set of indices $S\subseteq \{1,2,\ldots,n\}$ of maximum cardinality, such that the rows of $M_1$ (respectively $M_2$) indexed by $S$ are independent at the same time. This is a matroid intersection of two linear matroids defined by $M_1$ and $M_2$. We note that partition, graphic, and transversal matroids are special cases of linear matroids.

§ INDEPENDENCE-QUERY MATROID INTERSECTION ALGORITHM

In this section, we show that we can obtain an $\tO(nr^{3/4})$ matroid intersection algorithm in the dynamic-independence-oracle model.
This matches the state-of-the-art traditional independence-query algorithm of Blikstad [15].
We will only provide a proof sketch here because our algorithm is mostly an implementation of [15] in the new model with the help of (circuit) binary search trees.

Using the same construction as <ref> and <ref>, circuit binary search trees work verbatim in the dynamic-independence-oracle model (however, co-circuit binary search trees do not).
In particular, <ref><ref> can be checked with a single independence query.

For any integer $\beta \geq 1$, there exists a data structure that supports the following operations.

* $\textsc{Initialize}(\cM, S, Q_S, X)$: Given $S \in \cI$, the query-set $Q_S$ that corresponds to $S$, and $X \subseteq S$ or $X = \{t\}$, initialize the data structure in $\tO(|X|)$ time. The data structure also maintains $S$.

* $\textsc{Find}(y)$: Given $y \in \bar{S}$,

* if $X \subseteq S$, then return an $x \in X$ such that $S - y + x \in \cI$, or

* if $X = \{t\}$, then return the only element $x = s$ or $x = t$ in $X$ if $S + y \in \cI$ and $\bot$ otherwise.

The procedure returns $\bot$ if such an $x$ does not exist.
The procedure takes $\tO(\beta)$ time if the result is not $\bot$, and $\tO(1)$ time otherwise.

* $\textsc{Delete}(x)$: Given $x \in X$, if $x \not\in \{s, t\}$, delete $x$ from $X$ in $O(\log{n})$ time.

* $\textsc{Replace}(x, y)$: Given $x \in X$ and $y \not\in X$, replace $x$ in $X$ by $y$ in $O(\log{n})$ time.

* $\textsc{Update}(\Delta)$: Update $S$ to $S \oplus (\Delta \setminus \{s, t\})$ in amortized $\tO(\frac{|X| \cdot |\Delta|}{\beta})$ time.

The algorithm of [15] consists of the following three phases.

* First, obtain an $(1 - \epsilon)$-approximate solution $S$ using augmenting sets in $\tO(\frac{n\sqrt{r}}{\epsilon})$ time.

* Eliminate all augmenting paths in $G(S)$ of length at most $d$ using Cunningham's algorithm implemented by [24] in $\tO(nd + nr\epsilon)$ time.

* Finding the remaining $O(r/d)$ augmenting paths one at a time, using $\tO(n\sqrt{r})$ time each.

With $\epsilon = r^{-1/4}$ and $d = r^{3/4}$, the total running time is $\tO(nr^{3/4})$.
We briefly sketch how to implement the above three steps in the same running time also for the dynamic-independence-oracle model.

Note that the primary difficulty independence-query algorithms face is that we are only capable of checking <ref><ref> (using <ref>), which means that we can only explore the neighbors of $u \in \bar{S}$.
The aforementioned rank-query algorithms for building distance layers and finding blocking-flow style augmenting paths are thus inapplicable in the dynamic-independence-oracle model.

Approximation Algorithm.
The $O(n\sqrt{r}/\epsilon)$-query $(1-\epsilon)$-approximation algorithm of [15] needs to first compute distance layers up to distance $O(\frac{1}{\epsilon})$. This is done in a similar way as sketched below for “Eliminating Short Augmenting Paths”.

Otherwise, the approximation algorithm works through a series of “refine” operations
(algorithms , , and  in [24, 15]) to build a partial augmenting set. In each such operation, we only need to be able to do the following for some sets $(P,Q)$: start from some set $Q$ and find a maximal set $X\subseteq P$ such that
$Q+X$ is independent. This can be performed with a greedy algorithm in time (and dynamic query) $O(|P|)$, given that we already have queried set $Q$ before (which will be the case).

Finally, the approximation algorithm falls back to finding a special type of augmenting paths with respect to the current augmenting set, in the  algorithm of [15], with $\tO(n)$ queries for each such path. This algorithm can be implemented also in the dynamic-oracle model with the same query complexity.  relies on the  and  algorithms (which we already covered), in addition to a binary search trick to find feasible exchange pairs. This binary search can be implemented with the circuit trees (<ref>), and it takes a total time $\tO(n)$ to build them (since we can keep track of a queried set for $S$, and then we only need to build a circuit tree statically once for each layer in cost proportional to the size of the layer—which sums up to $\tO(n)$).

Eliminating Short Augmenting Paths.
Using [24]'s implementation of Cunningham's algorithm, we can eliminate all augmenting paths of length $d$, thereby obtaining a $(1 - 1/d)$-approximate solution.
The algorithm relies on <ref> to “fix” the distance layers after each augmentation.
Initially, all elements have distance $1$ or $2$ from $s$ depending on whether it belongs to $S$ (the common independent set obtained by the above approximation algorithm) or not.
Before the first and after each of the remaining $O(\epsilon r)$ augmentation, we can fix the distance layers as follows.

* For each $1 \leq \ell \leq d$ and $u \in L_{\ell}$, if $u$ is not of distance $\ell$ from $s$, i.e., there is no in-edge from $L_{\ell - 1}$ to $u$ anymore, move $u$ from $L_{\ell}$ to $L_{\ell + 2}$. This check is done as follows, depending on the parity of $\ell$.

* If $\ell$ is even, then for each $v \in L_{\ell - 1}$, we find all the unmarked $u \in L_{\ell}$ that $v$ has an edge to and mark $u$.
In the end, all the unmarked $u \in L_{\ell}$ do not belong to $L_{\ell}$ and should be moved to $L_{\ell + 2}$.

* If $\ell$ is odd, then we simply check if there is an in-edge from $L_{\ell - 1}$ to decide whether $u$ should be moved to other layers.

Both cases can be implemented efficiently with circuit binary search trees of <ref>: Each time we spend $\tO(1)$ time to either confirm that $u$ has distance $\ell$ from $s$ with respect to the current $S$ (in which case it will not be moved anymore in this iteration), or we increase the distance estimate of $u$.
The total running time is thus $\tO(nd + nr\epsilon)$, where $\tO(nd)$ comes from increasing the distance estimate of each element to at most $d$, and $\tO(nr\epsilon)$ comes from confirming that each element belongs to its distance layer in the $O(\epsilon r)$ iterations.

A caveat here is that we need to support insertion/deletion into the binary search trees.
This can be made efficient by doubling the size of a binary search tree (and re-initializing it) every time when there are not enough leaf nodes left in it.
The cost of re-building will be amortized to $\tO(1)$ time per update (i.e., movement of an element to another layer).

Finding a Single Augmenting Path.
With the $(1 - 1/d)$-approximate solution obtained in the first two steps, [15] then finds the remaining $O(r/d)$ augmenting paths one at a time, using the reachability algorithm of [17].
The reachability algorithm roughly goes as follows.
First, we initialize two circuit binary search trees (<ref>) over the two matroids for discovering out-edges and in-edges of elements in $\bar{S}$.
We then repeatedly run the following three steps until either an $(s, t)$-path is found (an arbitrary $(s, t)$-path suffices since such a path can be converted into a chordless one in $\tO(r)$ time along which augmentation is valid) or we conclude that $t$ is unreachable from $s$.
We keep track of a set of visited vertices $F$ which we know are reachable from $s$.

* Identify the set of unvisited heavy vertices in $\bar{S}$ that have at least $\sqrt{r}$ unvisited out-neighbors or has a direct edge toward $t$. This is done by sampling a set $R$ of unvisited vertices in $S$ and then computing for each vertex $u$ whether $R \cap \mathsf{OutNgh}(u) = \emptyset$, or equivalently, whether $S - R + u \in \cI$.
Intuitively, vertices with more out-neighbors are more likely to fail the test.
This can be tested for a single $R$ and all $u$ in $O(n)$ time in the dynamic-independence-oracle model.
With $O(\log{n})$ samples, heavy vertices can be successfully identified with high probability.

* Discover all the out-neighbors for each light vertex, taking a total of $\tO(n\sqrt{r})$ time using the circuit binary search tree over the whole run of the algorithm. (Each vertex turns from heavy to light at most once.)

* Perform a reversed breadth-first-search from all the heavy vertices simultaneously. We can assume that every vertex on the path is light (i.e., we find a “closest” heavy vertex reachable from $s$), and thus all its out-neighbors have already been discovered.
That is, going backward from $S$ to $\bar{S}$, we use the out-edges of light vertices.
From $\bar{S}$ to $S$, we use the circuit binary search tree.
This takes $\tO(n)$ time, and we either find a heavy vertex reachable from $s$ (in which case we make progress by visiting at least $\sqrt{r}$ vertices in $S$), or we conclude that all heavy vertices are unreachable from $s$ (in which case $t$ is unreachable either).

The number of iterations is bounded by $O(\sqrt{r})$ since we discover at least $\sqrt{r}$ unvisited vertices in $S$ every time.
The total running time of finding a single augmenting path is thus $\tO(n\sqrt{r})$.

Using the same parameters $\epsilon$ and $d$ as in [15] to combine the three phases, we obtain the following matroid intersection algorithm in the dynamic-independence-oracle model.

For two matroids $\cM_1 = (\ground, \cI_1)$ and $\cM_2 = (\ground, \cI_2)$, it takes $\tO(nr^{3/4})$ time to obtain the largest $S \in \cI_1 \cap \cI_2$ with high probability in the dynamic-independence-oracle model.