0.75pt 1.5pt

# Efficient Algorithms for Geometric Partial Matching††thanks: Work on this
paper was supported by NSF under grants CCF-15-13816, CCF-15-46392, and
IIS-14-08846, by an ARO grant W911NF-15-1-0408, and by BSF Grant 2012/229 from
the U.S.-Israel Binational Science Foundation.

Pankaj K. Agarwal  Department of Computer Science, Duke University. Email:
{pankaj, hc252<EMAIL_ADDRESS>Hsien-Chih Chang22footnotemark: 2 Allen
Xiao22footnotemark: 2

###### Abstract

Let $A$ and $B$ be two point sets in the plane of sizes $r$ and $n$
respectively (assume $r\leq n$), and let $k$ be a parameter. A matching
between $A$ and $B$ is a family of pairs in $A\times B$ so that any point of
$A\cup B$ appears in at most one pair. Given two positive integers $p$ and
$q$, we define the cost of matching $M$ to be $c(M)=\sum_{(a,b)\in
M}\mathopen{\|}a-b\mathclose{\|}_{p}^{q}$ where
$\mathopen{\|}\cdot\mathclose{\|}_{p}$ is the $L_{p}$-norm. The geometric
partial matching problem asks to find the minimum-cost size-$k$ matching
between $A$ and $B$.

We present efficient algorithms for geometric partial matching problem that
work for any powers of $L_{p}$-norm matching objective: An exact algorithm
that runs in $O((n+k^{2})\mathop{\mathrm{polylog}}n)$ time, and a
$(1+\varepsilon)$-approximation algorithm that runs in
$O((n+k\sqrt{k})\mathop{\mathrm{polylog}}n\cdot\log\varepsilon^{-1})$ time.
Both algorithms are based on the primal-dual flow augmentation scheme; the
main improvements involve using dynamic data structures to achieve efficient
flow augmentations. With similar techniques, we give an exact algorithm for
the planar transportation problem running in
$O(\min\\{n^{2},rn^{3/2}\\}\mathop{\mathrm{polylog}}n)$ time.

##  1 Introduction

Given two point sets $A$ and $B$ in the plane, we consider the problem of
finding the minimum-cost partial matching between $A$ and $B$. Formally,
suppose $A$ has size $r$ and $B$ has size $n$ where $r\leq n$. Let $G(A,B)$ be
the undirected complete bipartite graph between $A$ and $B$, and let the cost
of edge $(a,b)$ be
${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{$c(a,b)$}}}=\mathopen{\|}a-b\mathclose{\|}_{p}^{q}$,
for some positive integers $p$ and $q$. A _matching_ $M$ in $G(A,B)$ is a set
of edges sharing no endpoints. The _size_ of $M$ is the number of edges in
$M$. The cost of matching $M$, denoted _$c(M)$_ , is defined to be the sum of
costs of edges in $M$. For a parameter $k$, the problem of finding the
minimum-cost size-$k$ matching in $G(A,B)$ is called the _geometric partial
matching problem_. We call the corresponding problem in general bipartite
graphs (with arbitrary edge costs) the _partial matching problem_.111Partial
matching is also called _imperfect matching_ or _imperfect assignment_ [17,
9].

We also consider the following generalization of bipartite matching. Let
$\phi:A\cup B\to\mathbb{Z}$ be an integral _supply-demand function_ with
positive value on points of $A$ and negative value on points of $B$,
satisfying $\sum_{a\in A}\phi(a)=-\sum_{b\in B}\phi(b)$. Let
${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color<EMAIL_ADDRESS>A\cup B}\mathopen{|}\phi(p)\mathclose{|}$. A _transportation map_ is a
function $\tau:A\times B\to\mathbb{R}_{\geq 0}$ such that $\sum_{b\in
B}\tau(a,b)=\phi(a)$ for all $a\in A$ and $\sum_{a\in A}\tau(a,b)=-\phi(b)$
for all $b\in B$. We define the cost of $\tau$ to be

${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{$c(\tau)$}}}\coloneqq\sum_{(a,b)\in
A\times B}c(a,b)\cdot\tau(a,b).$

The _transportation problem_ asks to compute a transportation map of minimum
cost.

#### Related work.

Maximum-size bipartite matching is a classical problem in graph algorithms.
Upper bounds include the $O(m\sqrt{n})$ time algorithm by Hopcroft and Karp
[10] and the $O(m\min\\{\sqrt{m},n^{2/3}\\})$ time algorithm by Even and
Tarjan [7], where $n$ is the number of nodes and $m$ is the number of edges.
The first improvement in over thirty years was made by Mądry [15], which uses
an interior-point algorithm, runs in $O(m^{10/7}\mathop{\mathrm{polylog}}n)$
time.

The Hungarian algorithm [13] computes a minimum-cost maximum matching in a
bipartite graph in roughly $O(mn)$ time. Faster algorithms have been
developed, such as the $O(m\sqrt{n}\log(nC))$ time algorithms by Gabow and
Tarjan [8] and the improved $O(m\sqrt{n}\log C)$ time algorithm by Duan et al.
[6] assuming the edge costs are integral; here $C$ is the maximum cost of an
edge. Ramshaw and Tarjan [17] showed that the Hungarian algorithm can be
extended to compute a minimum-cost partial matching of size $k$ in
$O(km+k^{2}\log r)$ time, where $r$ is the size of the smaller side of the
bipartite graph. They also proposed a cost-scaling algorithm for partial
matching that runs in time $O(m\sqrt{k}\log(kC))$, again assuming that costs
are integral. By reduction to unit-capacity min-cost flow, Goldberg et al. [9]
developed a cost-scaling algorithm for partial matching with an identical
running time $O(m\sqrt{k}\log(kC))$, again only for integral edge costs.

In geometric settings, the Hungarian algorithm can be implemented to compute
an optimal perfect matching between $A$ and $B$ (assuming equal size) in time
$O(n^{2}\mathop{\mathrm{polylog}}n)$ [11] (see also [21, 1]). This algorithm
computes an optimal size-$k$ matching in time
$O(kn\mathop{\mathrm{polylog}}n)$. Faster approximation algorithms have been
developed for computing perfect matchings in geometric settings [21, 22, 4,
18]. Recall that the cost of the edges are the $q$th power of their
$L_{p}$-distances. When $q=1$, the best algorithm to date by Sharathkumar and
Agarwal [19] computes $(1+\varepsilon)$-approximation to the value of optimal
perfect matching in
$O(n\mathop{\mathrm{polylog}}n\cdot\mathop{\mathrm{poly}}\varepsilon^{-1})$
expected time with high probability. Their algorithm can also compute a
$(1+\varepsilon)$-approximate partial matching within the same time bound. For
$q>1$, the best known approximation algorithm to compute a perfect matching
runs in $O(n^{3/2}\mathop{\mathrm{polylog}}n\log(1/\varepsilon))$ time [18];
it is not obvious how to extend this algorithm to the partial matching
setting.

The transportation problem can also be formulated as an instance of the
minimum-cost flow problem. The strongly polynomial uncapacitated min-cost flow
algorithm by Orlin [16] solves the transportation problem in $O((m+n\log
n)n\log n)$ time. Lee and Sidford [14] give a weakly polynomial algorithm that
runs in $O(m\sqrt{n}\mathop{\mathrm{polylog}}(n,U))$ time, where $U$ is the
maximum amount of node supply-demand. Agarwal et al. [2, 3] showed that
Orlin’s algorithm can be implemented to solve 2D transportation in time
$O(n^{2}\mathop{\mathrm{polylog}}n)$. In case of $O(1)$-dimension Euclidean
space, by adapting the Lee-Sidford algorithm, they developed a
$(1+\varepsilon)$-approximation algorithm that runs in
$O(n^{3/2}\mathop{\mathrm{poly}}\varepsilon^{-1}\mathop{\mathrm{polylog}}(n,U))$
time. They also gave a Monte-Carlo algorithm that computes an
$O(\log^{2}(1/\varepsilon))$-approximate solution in $O(n^{1+\varepsilon})$
time with high probability. Recently, Khesin, Niklov, and Paramonov [12]
obtained a $(1+\varepsilon)$-approximation in low-dimensional Euclidean space
that runs in randomized
$O(n\mathop{\mathrm{poly}}\varepsilon^{-1}\mathop{\mathrm{polylog}}(n,U))$
time.

#### Our results.

There are three main results in this paper. First in Section 2 we present an
efficient algorithm for computing an optimal partial matching in the plane.

######  Theorem 1.1.

Given two point sets $A$ and $B$ in the plane each of size at most $n$ and an
integer $k\leq n$, a minimum-cost matching of size $k$ between $A$ and $B$ can
be computed in $O((n+k^{2})\mathop{\mathrm{polylog}}n)$ time.

We use _bichromatic closest pair (BCP)_ data structures to implement the
Hungarian algorithm efficiently, similar to Agarwal et al. [1] and Kaplan et
al. [11]. But unlike their algorithms which take $\Omega(n)$ time to find an
augmenting path, we show that after $O(n\mathop{\mathrm{polylog}}n)$ time
preprocessing, an augmenting path can be found in
$O(k\mathop{\mathrm{polylog}}n)$ time. The key is to recycle (rather than
rebuild) our data structures from one augmentation to the next. We refer to
this idea as the _rewinding mechanism_.

Next in Sections 3, we obtain a $(1+\varepsilon)$-approximation algorithm for
the geometric partial matching problem in the plane by providing an efficient
implementation of the unit-capacity min-cost flow algorithm by Goldberg et al.
[9].

######  Theorem 1.2.

Given two point sets $A$ and $B$ in $\mathbb{R}^{2}$ each of size at most $n$,
an integer $k\leq n$, and a parameter $\varepsilon>0$, a
$(1+\varepsilon)$-approximate min-cost matching of size $k$ between $A$ and
$B$ can be computed in
$O((n+k\sqrt{k})\mathop{\mathrm{polylog}}n\cdot\log\varepsilon^{-1})$ time.

The main challenge here is how to deal with the _dead nodes_ , which neither
have excess/deficit nor have flow passing through them, but still contribute
to the size of the graph. We show that the number of _alive nodes_ is only
$O(k)$, and then represent the dead nodes implicitly so that the Hungarian
search and computation of a blocking flow can be implemented in
$O(k\mathop{\mathrm{polylog}}n)$ time.

Finally in Section 4 we present a faster algorithm for the transportation
problem in $\mathbb{R}^{2}$ when the two point sets are unbalanced.

######  Theorem 1.3.

Given two point sets $A$ and $B$ in $\mathbb{R}^{2}$ of sizes $r$ and $n$
respectively with $r\leq n$, along with supply-demand function $\phi:A\cup
B\to\mathbb{Z}$, an optimal transportation map between $A$ and $B$ can be
computed in $O(\min\\{n^{2},rn^{3/2}\\}\mathop{\mathrm{polylog}}n)$ time.

Our result improves over the $O(n^{2}\mathop{\mathrm{polylog}}n)$ time
algorithm by Agarwal et al. [3] for $r=o(\sqrt{n})$. Similar to their
algorithm, we also use the strongly polynomial uncapacitated minimum-cost flow
algorithm by Orlin [16], but additional ideas are needed for efficient
implementation. Unlike in the case of matchings, the support of the
transportation problem may have size $\Omega(n)$ even when $r$ is a constant;
so naïvely we can no longer spend time proportional to the size of support of
the transportation map. However, with careful implementation we ensure that
the support is acyclic, and one can find an augmenting path in
$O(r\sqrt{n}\mathop{\mathrm{polylog}}n)$ time with proper data structures,
assuming $r\leq\sqrt{n}$.

##  2 Minimum-cost partial matchings using Hungarian algorithm

In this section, we solve the geometric partial matching problem and prove
Theorem 1.1 by implementing the Hungarian algorithm for partial matching in
$O((n+k^{2})\mathop{\mathrm{polylog}}n)$ time.

A node $v$ is _matched_ by matching $M$ if $v$ is the endpoint of some edge in
$M$; otherwise $v$ is _unmatched_. Given a matching $M$, an _augmenting path_
$\Pi=(a_{1},b_{1},\ldots,a_{\ell},b_{\ell})$ is an odd-length path with
unmatched endpoints ($a_{1}$ and $b_{\ell}$) that alternates between edges
outside and inside of $M$. The symmetric difference $M\oplus\Pi$ creates a new
matching of size $\mathopen{|}M\mathclose{|}+1$, called the _augmentation_ of
$M$ by $\Pi$. The dual to the standard linear program for partial matching has
one dual variable for each node $v$, called the _potential $\pi(v)$_ of $v$.
Given potential $\pi$, we can define the _reduced cost_ of the edges to be
${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{$c_{\pi}(v,w)$}}}\coloneqq
c(v,w)-\pi(v)+\pi(w)$. Potential $\pi$ is _feasible_ on edge $(v,w)$ if
$c_{\pi}(v,w)$ is nonnegative. Potential $\pi$ is _feasible_ if $\pi$ are
feasible on every edge in $G$. We say that an edge $(v,w)$ is _admissible_
under potential $\pi$ if $c_{\pi}(v,w)=0$.

#### Fast implementation of Hungarian search.

The Hungarian algorithm is initialized with $M\leftarrow\emptyset$ and
$\pi\leftarrow 0$. Each iteration of the Hungarian algorithm augments $M$ by
an admissible augmenting path $\Pi$, discovered using a procedure called the
_Hungarian search_. The algorithm terminates after $k$ augmentations, exactly
when $\mathopen{|}M\mathclose{|}=k$; Ramshaw and Tarjan [17] showed that $M$
is guaranteed to be an optimal partial matching.

The Hungarian search grows a set of _reachable nodes $X$_ from all unmatched
$v\in A$ using augmenting paths of admissible edges. Initially, $X$ is the set
of unmatched nodes in $A$. Let the _frontier_ of $X$ be the edges in $(A\cap
X)\times(B\setminus X)$. $X$ is grown by _relaxing_ an edge $(a,b)$ in the
frontier: Add $b$ into $X$, modify potential to make $(a,b)$ admissible,
preserve $c_{\pi}$ on other edges within $X$, and keep $\pi$ feasible on edges
outside of $X$. Specifically, the algorithm relaxes the min-reduced-cost
frontier edge $(a,b)$, and then raises $\pi(v)$ by $c_{\pi}(a,b)$ for all
$v\in X$. If $b$ is already matched, then we also relax the matching edge
$(a^{\prime},b)$ and add $a^{\prime}$ into $X$. The search finishes when $b$
is unmatched, and an admissible augmenting path now can be recovered.

In the geometric setting, we find the min-reduced-cost frontier edge using a
dynamic _bichromatic closest pair_ (BCP) data structure, similar to [3, 21].
Given two point sets $P$ and $Q$ in the plane and a weight function
$\omega:P\cup Q\to\mathbb{R}$, the BCP is two points $a\in P$ and $b\in Q$
minimizing the additively weighted distance $c(a,b)-\omega(a)+\omega(b)$.
Thus, a minimum reduced-cost frontier edge is precisely the BCP of point sets
$P=A\cap X$ and $Q=B\setminus X$, with $\omega=\pi$. Note that the “state” of
this BCP is parameterized by $X$ and $\pi$.

The dynamic BCP data structure by Kaplan et al. [11] supports point insertions
and deletions in $O(\mathop{\mathrm{polylog}}n)$ time and answers queries in
$O(\log^{2}n)$ time for our setting. Each relaxation in the Hungarian search
requires one query, one deletion, and at most one insertion (aside from the
potential updates). As $\mathopen{|}M\mathclose{|}\leq k$ throughout, there
are at most $2k$ relaxations in each Hungarian search, and the BCP can be used
to implement each Hungarian search in $O(k\mathop{\mathrm{polylog}}n)$ time.

#### Rewinding mechanism.

We observe that exactly one node of $A$ is newly matched after an
augmentation. Thus (modulo potential changes), we can obtain the initial state
of the BCP for the $(i+1)$-th Hungarian search from the $i$-th one with a
single BCP deletion.

If we remember the sequence of points added to $X$ in the $i$-th Hungarian
search, then at the start of the $(i+1)$-th Hungarian search we can _rewind_
this sequence by applying the opposite insert/delete operation to each BCP
update in reverse order to obtain the initial state of the $i$-th BCP. With
one additional BCP deletion, we have the initial state of the $(i+1)$-th BCP.
The number of insertions/deletions is bounded by the number of relaxations per
Hungarian search which is $O(k)$. Therefore we can recover, in
$O(k\mathop{\mathrm{polylog}}n)$ time, the initial BCP data structure for each
Hungarian search beyond the first. We refer to this procedure as the
_rewinding mechanism_.

#### Potential updates.

We modify a trick from Vaidya [21] to batch potential updates. Potential is
tracked with a _stored value_ $\gamma(v)$, while the _true value_ of $\pi(v)$
may have changed since $\gamma(v)$ was last recorded. This is done by
aggregating potential changes into a variable $\delta$, which is initially 0
at the very beginning of the algorithm. Whenever we would raise the potential
of all nodes in $X$, we raise $\delta$ by that amount instead. We maintain the
following invariant: $\pi(v)=\gamma(v)$ for $v\not\in X$, and
$\pi(v)=\gamma(v)+\delta$ for $v\in X$.

At the beginning of the algorithm, $X$ is empty and stored values are equal to
true values. When $a\in A$ is added to $X$, we update its stored value to
$\pi(a)-\delta$ for the current value of $\delta$, and use that stored value
as its BCP weight. Since the BCP weights are uniformly offset from $\pi(v)$ by
$\delta$, the pair reported by the BCP is still minimum. When $b\in B$ is
added to $X$, we update its stored value to $\pi(b)-\delta$ (although it won’t
be added to a BCP set). When a node is removed from $X$ (e.g. by augmentation
or rewinding), we update the stored potential
$\gamma(v)\leftarrow\pi(v)+\delta$, again for the current value of $\delta$.
Unlike Vaidya [21], we do not reset $\delta$ across Hungarian searches for the
sake of rewinding.

There are $O(k)$ relaxations and thus $O(k)$ updates to $\delta$ per Hungarian
search. $O(k)$ stored values are updated per rewinding, so the time spent on
potential updates per Hungarian search is $O(k)$. Putting everything together,
our implementation of the Hungarian algorithm runs in
$O((n+k^{2})\mathop{\mathrm{polylog}}n)$ time. This proves Theorem 1.1.

##  3 Approximating min-cost partial matching through cost-scaling

In this section we describe an approximation algorithm for computing a min-
cost partial matching. We reduce the problem to computing a min-cost
circulation in a flow network (Section 3.1). We adapt the cost-scaling
algorithm by Goldberg et al. [9] for computing min-cost flow of a unit-
capacity network (Section 3.2). Finally, we show how their algorithm can be
implemented in
$O\left((n+k^{3/2})\mathop{\mathrm{polylog}}(n)\log(1/\varepsilon)\right)$
time in our setting (Section 3.3).

###  3.1 From matching to circulation

Given a bipartite graph $G$ with node sets $A$ and $B$, we construct a flow
network $N=(V,\vec{E})$ in a standard way [17] so that a min-cost matching in
$G$ corresponds to a min-cost integral circulation in $N$.

#### Flow network.

Each node in $G$ becomes a node in $N$ and each edge $(a,b)$ in $G$ becomes an
arc $a\mathord{\shortrightarrow}b$ in $N$; we refer to these nodes and arcs as
_bipartite nodes_ and _bipartite arcs_. We also include a _source_ node $s$
and _sink_ node $t$ in $N$. For each $a\in A$, we add a _left dummy arc_
$s\mathord{\shortrightarrow}a$ and for each $b\in B$ a _right dummy arc_
$b\mathord{\shortrightarrow}t$. The cost _$c(v\mathord{\shortrightarrow}w)$_
is equal to $c(v,w)$ if $v\mathord{\shortrightarrow}w$ is a bipartite arc and
$0$ if $v\mathord{\shortrightarrow}w$ is a dummy arc. All arcs in $N$ have
unit capacity.

Let $\phi:V\to\mathbb{Z}$ be an integral supply/demand function on nodes of
$N$. The positive values of $\phi(v)$ are referred to as _supply_ , and the
negative values of $\phi(v)$ as _demand_. A _pseudoflow_ $f:\vec{E}\to[0,1]$
is a function on arcs of $N$. The _support_ of $f$ in $N$, denoted as
_$\operatorname{supp}(f)$_ , is the set of arcs with positive flow. Given a
pseudoflow $f$, the _imbalance_ of a node is

${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{$\phi_{f}(v)$}}}\coloneqq\phi(v)+\sum_{w\mathord{\shortrightarrow}v\in\vec{E}}{f(w\mathord{\shortrightarrow}v)}-\sum_{v\mathord{\shortrightarrow}w\in\vec{E}}{f(v\mathord{\shortrightarrow}w)}.$

We call positive imbalance _excess_ and negative imbalance _deficit_. A node
is _balanced_ if it has zero imbalance. If all nodes are balanced, the
pseudoflow is a _circulation_. The _cost_ of a pseudoflow is defined to be

${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{$c(f)$}}}\coloneqq\sum_{v\mathord{\shortrightarrow}w\in\operatorname{supp}(f)}c(v\mathord{\shortrightarrow}w)\cdot
f(v\mathord{\shortrightarrow}w).$

The _minimum-cost flow problem_ (MCF) asks to find a circulation of minimum
cost.

If we set $\phi(s)=k$, $\phi(t)=k$, and $\phi(v)=0$ for all $v\in A\cup B$,
then an integral circulation $f$ corresponds to a partial matching $M$ of size
$k$ and vice versa. Moreover, $c(M)=c(f)$. Hence, the problem of computing a
min-cost matching of size $k$ in $G$ transforms to computing an integral
circulation in $N$. The following lemma will be useful for our algorithm.

######  Lemma 3.1.

Let $N$ be the network constructed from the bipartite graph $G$ above.

  1. (i)

For any integral circulation $g$ in $N$, the size of $\operatorname{supp}(g)$
is at most $3k$.

  2. (ii)

For any integral pseudoflow $f$ in $N$ with $O(k)$ excess, the size of
$\operatorname{supp}(f)$ is $O(k)$.

###  3.2 A cost-scaling algorithm

Before describing the algorithm, we need to introduce a few more concepts.

#### Residual network and admissibility.

If $f$ is an integral pseudoflow on $N$ (that is,
$f(v\mathord{\shortrightarrow}w)\in\\{0,1\\}$ for every arc in $\vec{E}$),
then each arc $v\mathord{\shortrightarrow}w$ in $N$ is either _idle_ with
$f(v\mathord{\shortrightarrow}w)=0$ or _saturated_ with
$f(u\mathord{\shortrightarrow}v)=1$.

Given a pseudoflow $f$, the _residual network_ $N_{f}=(V,\vec{E}_{f})$ is
defined as follows. For each idle arc $v\mathord{\shortrightarrow}w$ in
$\vec{E}$, we add a _forward_ residual arc $v\mathord{\shortrightarrow}w$ in
$N_{f}$. For each saturated arc $v\mathord{\shortrightarrow}w$ in $\vec{E}$,
we add a _backward_ residual arc $w\mathord{\shortrightarrow}v$ in $N_{f}$.
The set of residual arcs in $N_{f}$ is therefore

$\vec{E}_{f}\coloneqq\\{v\mathord{\shortrightarrow}w\mid
f(v\mathord{\shortrightarrow}w)=0\\}\cup\\{w\mathord{\shortrightarrow}v\mid
f(v\mathord{\shortrightarrow}w)=1\\}.$

The cost of a forward residual arc $v\mathord{\shortrightarrow}w$ is
$c(v\mathord{\shortrightarrow}w)$, while the cost of a backward residual arc
$w\mathord{\shortrightarrow}v$ is $-c(v\mathord{\shortrightarrow}w)$. Each arc
in $N_{f}$ also has unit capacity. By Lemma 3.1, $N_{f}$ has $O(k)$ backward
arcs if $f$ has $O(k)$ excess.

A _residual pseudoflow_ $g$ in $N_{f}$ can be used to change $f$ into a
different pseudoflow on $N$ by _augmentation_. For simplicity, we only
describe augmentation for the case where $f$ and $g$ are integral.
Specifically, augmenting $f$ by $g$ produces a pseudoflow $f^{\prime}$ in $N$
where

$f^{\prime}(v\mathord{\shortrightarrow}w)=\begin{cases}0&{w\mathord{\shortrightarrow}v}\in\vec{E}_{f}\text{
and }g(w\mathord{\shortrightarrow}v)=1,\\\
1&{v\mathord{\shortrightarrow}w}\in\vec{E}_{f}\text{ and
}g(v\mathord{\shortrightarrow}w)=1,\\\
f(v\mathord{\shortrightarrow}w)&\text{otherwise.}\end{cases}$

Using LP duality for min-cost flow, we assign _potential $\pi(v)$_ to each
node $v$ in $N$. The _reduced cost_ of an arc $v\mathord{\shortrightarrow}w$
in $N$ with respect to $\pi$ is defined as

$c_{\pi}(v\mathord{\shortrightarrow}w)\coloneqq
c(v\mathord{\shortrightarrow}w)-\pi(v)+\pi(w).$

Similarly we define the reduced cost of arcs in $N_{f}$: the reduced cost of a
forward residual arc $v\mathord{\shortrightarrow}w$ in $N_{f}$ is
$c_{\pi}(v\mathord{\shortrightarrow}w)$, and the reduced cost of a backward
residual arc $w\mathord{\shortrightarrow}v$ in $N_{f}$ is
$-c_{\pi}(v\mathord{\shortrightarrow}w)$. Abusing the notation, we also use
$c_{\pi}$ to denote the reduced cost of arcs in $N_{f}$.

The _dual feasibility constraint_ asks that
$c_{\pi}(v\mathord{\shortrightarrow}w)\geq 0$ holds for every arc
$v\mathord{\shortrightarrow}w$ in $\vec{E}$; potential $\pi$ that satisfy this
constraint is said to be _feasible_. Suppose we relax the dual feasibility
constraint to allow some small violation in the value of
$c_{\pi}(v\mathord{\shortrightarrow}w)$. We say that a pair of pseudoflow $f$
and potential $\pi$ is _$\theta$ -optimal_ [20, 5] if
$c_{\pi}(v\mathord{\shortrightarrow}w)\geq-\theta$ for every residual arc
$v\mathord{\shortrightarrow}w$ in $\vec{E}_{f}$. Pseudoflow $f$ is _$\theta$
-optimal_ if it is $\theta$-optimal with respect to some potential $\pi$;
potential $\pi$ is _$\theta$ -optimal_ if it is $\theta$-optimal with respect
to some pseudoflow $f$. Given a pseudoflow $f$ and potential $\pi$, a residual
arc $v\mathord{\shortrightarrow}w$ in $\vec{E}_{f}$ is _admissible_ if
$c_{\pi}(v\mathord{\shortrightarrow}w)\leq 0$. We say that a pseudoflow $g$ in
$N_{f}$ is _admissible_ if $g(v\mathord{\shortrightarrow}w)>0$ only on
admissible arcs $v\mathord{\shortrightarrow}w$, and
$g(v\mathord{\shortrightarrow}w)=0$ otherwise.222The same
admissibility/feasibility definitions will be used later in Section 4.
However, the algorithm in Section 4 maintains a 0-optimal $f$ and therefore
admissible residual arcs always have
$c_{\pi}(v\mathord{\shortrightarrow}w)=0$. We will use the following well-
known property of $\theta$-optimality.

######  Lemma 3.2.

Let $f$ be an $\theta$-optimal pseudoflow in $N$ and let $g$ be an admissible
pseudoflow in $N_{f}$. Then $f$ augmented by $g$ is also $\theta$-optimal in
$N$.

Using Lemma 3.1, the following lemma can be proved about $\theta$-optimality:

######  Lemma 3.3.

Let $f$ be a $\theta$-optimal integer circulation in $N$, and $f^{*}$ be an
optimal integer circulation for $N$. Then, $c(f)\leq c(f^{*})+6k\theta$.

#### Estimating the value of $c(f^{*})$.

We now describe a procedure for estimating $c(f^{*})$ within a polynomial
factor, which is used to initialize the cost-scaling algorithm.

Let _$T$_ be a minimum spanning tree of $A\cup B$ under the cost function $c$.
Let $e_{1},e_{2},\ldots,e_{n-1}$ be the edges of $T$ sorted in nondecreasing
order of length. Let _$T_{i}$_ be the forest consisting of the nodes of $A\cup
B$ and edges $e_{1},\ldots,e_{i}$. We call a matching $M$ _intra-cluster_ if
both endpoints of each edge in $M$ lie in the same connected component of
$T_{i}$. The following lemma will be used by our cost-scaling algorithm:

######  Lemma 3.4.

Let $A$ and $B$ be two point sets in the plane. Define _$i^{*}$_ to be the
smallest index $i$ such that there is an intra-cluster matching of size $k$ in
$T_{i^{*}}$. Set _$\overline{\theta}$_ $\coloneqq n^{q}\cdot c(e_{i^{*}})$.
Then

  1. (i)

The value of $i^{*}$ can be computed in $O(n\log n)$ time.

  2. (ii)

$c(e_{i^{*}})\leq c(f^{*})\leq\overline{\theta}$.

  3. (iii)

There is a $\overline{\theta}$-optimal circulation in the network $N$ with
respect to the all-zero potential, assuming $\phi(s)=k$, $\phi(t)=-k$, and
$\phi(v)=0$ for all $v\in A\cup B$.

As a consequence of Lemmas 3.4(ii) and 3.3, we have:

######  Corollary 3.5.

The cost of a $\underline{\theta}$-optimal integral circulation in $N$ is at
most $(1+\varepsilon)c(f^{*})$, where
${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color<EMAIL_ADDRESS>c(e_{i^{*}})$.

#### Overview of the algorithm.

We are now ready to describe our algorithm, which closely follows Goldberg et
al. [9]. The algorithm works in scales. In the beginning of each scale, we fix
a _cost scaling parameter_ _$\theta$_ and maintain potential $\pi$ with the
following property:

($\ast$)

There exists a $2\theta$-optimal integral circulation in $N$ with respect to
$\pi$.

For the initial scale, we set $\theta\leftarrow\overline{\theta}$ and
$\pi\leftarrow 0$. By Lemma 3.4(iii), property (($\ast$) ‣ 3.2) is satisfied
initially. Each scale of the algorithm consists of two stages. In the _scale
initialization_ stage, Scale-Init computes a $\theta$-optimal pseudoflow $f$.
In the _refinement_ stage, Refine converts $f$ into a $\theta$-optimal
(integral) circulation $g$. In both stages, $\pi$ is updated as necessary. If
$\theta\leq\underline{\theta}$, we return $g$. Otherwise, we set
$\theta\leftarrow\theta/2$ and start the next scale. Note that property
(($\ast$) ‣ 3.2) is satisfied in the beginning of each scale.

By Corollary 3.5, when the algorithm terminates, it returns an integral
circulation $\tilde{f}$ in $N$ of cost at most $(1+\varepsilon)c(f^{*})$,
which corresponds to a $(1+\varepsilon)$-approximate min-cost matching of size
$k$ in $G$. The algorithm terminates in
$\log_{2}(\overline{\theta}/\underline{\theta})=O(\log(n/\varepsilon))$
scales.

#### Scale initialization.

In the first scale, we compute a $\overline{\theta}$-optimal pseudoflow by
simply setting $f(v\mathord{\shortrightarrow}w)\leftarrow 0$ for all arcs in
$\vec{E}$. For subsequent scales, we begin with a $2\theta$-optimal
circulation $f$ in $N$. First, we raise the potential of all nodes in $A$ by
$\theta$, all nodes in $B$ by $2\theta$, and $t$ by $3\theta$. The potential
of $s$ is unchanged. Such potential change increases the reduced cost of all
forward arcs to at least $-\theta$.

Next, for each backward arc $w\mathord{\shortrightarrow}v$ in $N_{f}$ with
$c_{\pi}(w\mathord{\shortrightarrow}v)<-\theta$, we set
$f(v\mathord{\shortrightarrow}w)\leftarrow 0$ (that is, make arc
$v\mathord{\shortrightarrow}w$ idle), which replaces the backward arc
$w\mathord{\shortrightarrow}v$ in $N_{f}$ with forward arc
$v\mathord{\shortrightarrow}w$ of positive reduced cost. After this step, the
resulting pseudoflow must be $\theta$-optimal as all arcs of $N_{f}$ have
reduced cost at least $-\theta$.

The desaturation of each backward arc creates one unit of excess. Since there
are at most $3k$ backward arcs, the pseudoflow has at most $3k$ excess after
Scale-Init. There are $O(n)$ potential updates and $O(k)$ arcs to desautrate,
so the time required for Scale-Init is $O(n)$.

#### Refinement.

The procedure Refine converts a $\theta$-optimal pseudoflow with $O(k)$ excess
into a $\theta$-optimal circulation, using a primal-dual augmentation
algorithm. A path in $N_{f}$ is an _augmenting path_ if it begins at an excess
node and ends at a deficit node. We call an admissible pseudoflow $g$ in
$N_{f}$ an _admissible blocking flow_ if $g$ saturates at least one arc in
every admissible augmenting path in $N_{g}$. In other words, there is no
admissible excess-deficit path in the residual network after augmentation by
$g$. Each iteration of Refine finds an admissible blocking flow to be added to
the current pseudoflow in two steps:

  1. 1.

_Hungarian search_ : a Dijkstra-like search that begins at the set of excess
nodes and raises potential until there is an excess-deficit path of admissible
arcs in $N_{f}$.

  2. 2.

_Augmentation_ : construct an admissible blocking flow by performing depth-
first search on the set of admissible arcs of $N_{f}$. It suffices to
repeatedly extract admissible augmenting paths until no more admissible
excess-deficit paths remain.

The algorithm repeats these steps until the total excess becomes zero. The
following lemma bounds the number of iterations in the Refine procedure at
each scale.

######  Lemma 3.6.

Let $\theta$ be the scaling parameter and $\pi_{0}$ the potential function at
the beginning of a scale, such that there exists an integral $2\theta$-optimal
circulation with respect to $\pi_{0}$. Let $f$ be a $\theta$-optimal
pseudoflow with excess $O(k)$. Then Refine terminates within $O(\sqrt{k})$
iterations.

###### Proof.

We sketch the proof, which is adapted from Goldberg et al. [9]. Let $f_{0}$ be
the assumed $2\theta$-optimal integral circulation with respect to $\pi_{0}$,
and let $\pi$ be the potential maintained during Refine. Let
${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{$d(v)$}}}\coloneqq(\pi(v)-\pi_{0}(v))/\theta$,
that is, the increase in potential at $v$ in units of $\theta$. We divide the
iterations of Refine into two phases: before and after every (remaining)
excess node has $d(v)\geq\sqrt{k}$. Each Hungarian search raises excess
potential by at least $\theta$, since we use blocking flows. Thus, the first
phase lasts at most $\sqrt{k}$ iterations.

At the start of the second phase, consider the set of arcs
${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color<EMAIL_ADDRESS>f(v\mathord{\shortrightarrow}w)<f_{0}(v\mathord{\shortrightarrow}w)\\}$. One
can argue that the remaining excess with respect to $f$ is bounded above by
the size of any cut separating the excess and deficit nodes [9, Lemma 4]. The
proof examines cuts
${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color<EMAIL_ADDRESS>d(v)>i\\}$ for $0\leq i\leq\sqrt{k}$. By $\theta$-optimality of $f$ and
$2\theta$-optimality of $f_{0}$, one can show that each arc in $E^{+}$ crosses
at most 3 cuts. Furthermore, the size of $E^{+}$ is $O(k)$, bounded by the
support size of $f$ and $f_{0}$. Averaging, there is a cut among $Y_{i}$s of
size at most $3k/\sqrt{k}$, so the total excess remaining is $O(\sqrt{k})$.
Each iteration of Refine eliminates at least one unit of excess, so the number
of second phase iterations is also at most $O(\sqrt{k})$. ∎

In the next subsection we show that after $O(n\mathop{\mathrm{polylog}}n)$
time preprocessing, an iteration of Refine can be performed in
$O(k\mathop{\mathrm{polylog}}n)$ time (Lemma 3.8). By Lemma 3.6 and the fact
the algorithm terminates in $O(\log(n/\varepsilon))$ scales, the overall
running time of the algorithm is
$O((n+k^{3/2})\mathop{\mathrm{polylog}}n\log(1/\varepsilon))$, as claimed in
Theorem 1.2.

###  3.3 Fast implementation of refinement stage

We now describe a fast implementation of Refine. The Hungarian search and
augmentation steps are similar: each traversing through the residual network
using admissible arcs starting from the excess nodes. Due to lack of space, we
only describe the former.

At a high level, let _$X$_ be the subset of nodes visited by the Hungarian
search so far. Initially $X$ is the set of excess nodes. At each step, the
algorithm finds a minimum-reduced-cost arc $v\mathord{\shortrightarrow}w$ in
$N_{f}$ from $X$ to $V\setminus X$. If $v\mathord{\shortrightarrow}w$ is not
admissible, the potential of all nodes in $X$ is increased by $\left\lceil
c_{\pi}(v\mathord{\shortrightarrow}w)/\theta\right\rceil$ to make
$v\mathord{\shortrightarrow}w$ admissible. If $w$ is a deficit node, the
search terminates. Otherwise, $w$ is added to $X$ and the search continues.

Implementing the Hungarian search efficiently is more difficult than in
Section 2 because (a) excess nodes may show up in $A$ as well as in $B$, (b) a
balanced node may become imbalanced later in the scales, and (c) the potential
of excess nodes may be non-uniform. We therefore need a more complex data
structure.

We call a node $v$ of $N$ _dead_ if $\phi_{f}(v)=0$ and no arc of
$\operatorname{supp}(f)$ is incident to $v$; otherwise $v$ is _alive_. Note
that $s$ and $t$ are always alive. Let _${A}^{\text{\char 96}}$_ denote the
set of alive nodes in $A$; define _${B}^{\text{\char 96}}$_ similarly. There
are only $O(k)$ alive nodes, as each can be charged to its adjoining
$\operatorname{supp}(f)$ arcs or its imbalance. We treat alive and dead nodes
separately to implement the Hungarian search efficiently. By definition, dead
nodes only adjoin forward arcs in $N_{f}$. Thus, the in-degree (resp. out-
degree) of a node in ${A}\setminus{A}^{\text{\char 96}}$ (resp.
${B}\setminus{B}^{\text{\char 96}}$) is 1, and any path passing through a dead
node has a subpath of the form
$s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}b$ for some $b\in B$ or
$a\mathord{\shortrightarrow}v\mathord{\shortrightarrow}t$ for some $a\in A$.
Consequently, a path in $N_{f}$ may have at most two consecutive dead nodes,
and in the case of two consecutive dead nodes there is a subpath of the form
$s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}w\mathord{\shortrightarrow}t$
where $v\in{A}\setminus{A}^{\text{\char 96}}$ and
$w\in{B}\setminus{B}^{\text{\char 96}}$. We call such paths, from an alive
node to an alive node with only dead interior nodes, _alive paths_. Let the
reduced cost $c_{\pi}(\Pi)$ of an alive path $\Pi$ be the sum of $c_{\pi}$
over its arcs. We say $\Pi$ is _weakly admissible_ if $c_{\pi}(\Pi)\leq 0$.

We find the min-reduced-cost alive path of lengths $1$, $2$, and $3$ leaving
$X$, then relax the cheapest among them (raise potential of $X$ by $c(\Pi)$
and add every node of $\Pi$ into $X$). Essentially, relaxing alive paths
“skips over” dead nodes. Since reduced costs telescope on paths, weak
admissibility of an alive path depends only on the potential of its alive
endpoints. Thus, we can query the minimum alive path using a partial
assignment of $\pi$ on only the alive nodes, leaving $\pi$ over the dead nodes
untracked. We now describe a data structure for each path length. Note that
our “time budget” per Hungarian search is $O(k\mathop{\mathrm{polylog}}n)$.

#### Finding length-1 paths.

This data structure finds a min-reduced-cost arc from an alive node of $X$ to
an alive node of $V\setminus X$. There are $O(k)$ backward arcs, so the
minimum among backward arcs can be maintained explicitly in a priority queue
and retrieved in $O(1)$ time.

There are three types of forward arcs: $s\mathord{\shortrightarrow}a$ for some
$a\in{A}^{\text{\char 96}}$, $b\mathord{\shortrightarrow}t$ for some
$b\in{B}^{\text{\char 96}}$, and bipartite arc $a\mathord{\shortrightarrow}b$
with two alive endpoints. Arcs of the first (resp. second) type can be found
by maintaining ${A}^{\text{\char 96}}\setminus X$ (resp. ${B}^{\text{\char
96}}\cap X$) in a priority queue, but should only be queried if $s\in X$
(resp. $t\not\in X$). The cheapest arc of the third type can be maintained
using a dynamic (additively weighted) bichromatic closest pair (BCP) data
structure between ${A}^{\text{\char 96}}\cap X$ and ${B}^{\text{\char
96}}\setminus X$, with reduced cost as the weighted pair distance. Such BCP
data structure can be implemented so that insertions/deletions can be
performed in $O(\mathop{\mathrm{polylog}}k)$ time [11].

#### Finding length-2 paths.

We describe how to find a cheapest path of the form
$s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}b$ where $v$ is dead
and $b\in{B}^{\text{\char 96}}$. A cheapest path
$a\mathord{\shortrightarrow}v\mathord{\shortrightarrow}t$ can be found
similarly. Similar to length-1 paths, we only query paths starting at $s$ if
$s\in X$ and paths ending at $t$ if $t\not\in X$.

Note that
$c_{\pi}(s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}b)=c(v,b)+\pi(b)-\pi(s)$.
Since $\pi(s)$ is common in all such paths, it suffices to find a pair $(v,w)$
between ${A}\setminus{A}^{\text{\char 96}}$ and ${B}^{\text{\char
96}}\setminus X$ minimizing $c(v,w)+\pi(w)$. This is done by maintaining a
dynamic BCP data structure between ${A}\setminus{A}^{\text{\char 96}}$ and
${B}^{\text{\char 96}}\setminus X$ with the cost of a pair $(v,w)$ being
$c(v,w)+\pi(w)$. We may require an update operation for each alive node added
to $X$ during the Hungarian search, of which there are $O(k)$, so the time
spent during a search is $O(k\mathop{\mathrm{polylog}}n)$.

Since the size of ${A}\setminus{A}^{\text{\char 96}}$ is at least $r-k$, we
cannot construct this BCP from scratch at the beginning of each iteration. To
resolve this, we use the idea of rewinding from Section 2, with a slight
twist. There are now _two_ ways that the initial BCP may change across
consecutive Hungarian searches: (1) the initial set $X$ may change as nodes
lose excess through augmentation, and (2) the set of alive/dead nodes in $A$
may change. The first is identical to the situation in Section 2; the number
of excess depletions is $O(k)$ over the course of Refine. For the second, the
alive/dead status of a node can change only if the blocking flow found passes
through the node. By Lemma 3.7 below, there are $O(k)$ such changes per
Hungarian search, which can be done in $O(k\mathop{\mathrm{polylog}}n)$ time.

#### Finding length-3 paths.

We now describe how to find the cheapest path of the form
$s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}w\mathord{\shortrightarrow}t$
where $v\in{A}\setminus{A}^{\text{\char 96}}$ and
$w\in{B}\setminus{B}^{\text{\char 96}}$. Note that
$c_{\pi}(s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}w\mathord{\shortrightarrow}t)=c(v\mathord{\shortrightarrow}w)-\pi(s)+\pi(t)$.
A pair $(v,w)$ between ${A}\setminus{A}^{\text{\char 96}}$ and
${B}\setminus{B}^{\text{\char 96}}$ minimizing $c(v,w)$ can be found by
maintaining a dynamic BCP data structure similar to the case of length-2
paths.

This BCP data structure has no dependency on $X$—the only update required
comes from membership changes to ${A}^{\text{\char 96}}$ or ${B}^{\text{\char
96}}$ after an augmentation. Applying Lemma 3.7 again, there are $O(k)$
alive/dead updates caused by an augmentation, so the time for these updates
per Hungarian search is $O(k\mathop{\mathrm{polylog}}n)$.

#### Updating potential.

Potential updates for alive nodes can be handled in a batched fashion as in
Section 2. The three data structures above have no dependency on the dead node
potential; we leave them untracked as described before. The Hungarian search
remains intact since alive nodes are visited in the same order as when using
arc-by-arc relaxations. However, we need values of $\pi$ on all nodes at the
end of a scale (for the next Scale-Init) and for individual dead nodes
whenever they become alive (after augmentation).

We can reconstruct a “valid” potential in these situations. To recover
potential for $v\in{A}\setminus{A}^{\text{\char 96}}$ we set
$\pi(v)\leftarrow\pi(s)$, and for $v\in{B}\setminus{B}^{\text{\char 96}}$ we
set $\pi(v)\leftarrow\pi(t)$. Straightforward calculation shows that such
potential (1) preserves $\theta$-optimality, and (2) makes $\Pi$ (arc-wise)
admissible for any weakly admissible alive path $\Pi$. Hence, a blocking flow
composed of weakly admissible alive paths is admissible under the recovered
potential.

The following lemma is crucial to the analysis of running time for the
Hungarian search, bounding both the number of relaxations and potential
update/recovery operations.

######  Lemma 3.7.

Both Hungarian search and augmentation stages explore $O(k)$ nodes, and the
blocking flow found in augmentation stage is incident to $O(k)$ nodes.

Augmentation can also be implemented in $O(k\mathop{\mathrm{polylog}}n)$ time,
after $O(n\mathop{\mathrm{polylog}}n)$ time preprocessing, using similar data
structures. We thus obtain the following:

######  Lemma 3.8.

After $O(n\mathop{\mathrm{polylog}}n)$ time preprocessing, each iteration of
Refine can be implemented in $O(k\mathop{\mathrm{polylog}}n)$ time.

##  4 Transportation algorithm

Given two point sets $A$ and $B$ in $\mathbb{R}^{2}$ of sizes $r$ and $n$
respectively and a supply-demand function $\phi:A\cup B\to\mathbb{Z}$ as
defined in the introduction, we present an
$O(rn^{3/2}\mathop{\mathrm{polylog}}n)$ time algorithm for computing an
optimal transport map between $A$ and $B$. By applying this algorithm in the
case of $r\leq\sqrt{n}$ and the one by Agarwal et al. [3] when $r>\sqrt{n}$,
we prove Theorem 1.3. We use a standard reduction to the uncapacitated min-
cost flow problem and use Orlin’s algorithm [16] as well as some of the ideas
from Agarwal et al. [3] for efficient implementation under the geometric
settings. We first present an overview of the algorithm and then describe its
fast implementation that achieves the desired running time.

###  4.1 Overview of the algorithm

Orlin’s algorithm follows an excess-scaling paradigm and the primal-dual
framework. It maintains a _scale parameter_ $\Delta$, a flow function $f$, and
potential $\pi$ on the nodes. Initially $\Delta$ is equal to the total supply,
$f=0$, and $\pi=0$. We fix a constant parameter $\alpha\in(0.5,1)$. A node $v$
is called _active_ if the magnitude of imbalance of $v$ is at least
$\alpha\Delta$. At each step, using the Hungarian search, the algorithm finds
an admissible excess-to-deficit path between active nodes in the residual
network and pushes a flow of amount $\Delta$ along this path.333Note that this
augmentation may convert an excess node into a deficit node. Repeat the
process until either active excess or deficit nodes are gone; when this
happens, $\Delta$ is halved. The sequence of augmentations with a fixed value
of $\Delta$ is called an _excess scale_.

The algorithm also performs two preprocessing steps at the beginning of each
excess scale. First, if $f(v\mathord{\shortrightarrow}w)\geq 3n\Delta$,
$v\mathord{\shortrightarrow}w$ is contracted to a single node $z$ with
$\phi(z)=\phi(v)+\phi(w)$.444Intuitively, $f(v\mathord{\shortrightarrow}w)$ is
so high that future scales cannot deplete the flow on
$v\mathord{\shortrightarrow}w$. Second, if there are no active excess nodes
and $f(v\mathord{\shortrightarrow}w)=0$ for every arc
$v\mathord{\shortrightarrow}w$, then $\Delta$ is aggresively lowered to
$\max_{v}\phi(v)$.

When the algorithm terminates, an optimal circulation in the contracted
network is found. We use the algorithm described in Agarwal et al. [3] to
recover an optimal circulation for the original network in
$O(n\mathop{\mathrm{polylog}}n)$ time. Orlin showed that the algorithm
terminates within $O(n\log n)$ scales and performs a total of $O(n\log n)$
augmentations. In the next subsection, we describe an algorithm that, after
$O(n\mathop{\mathrm{polylog}}n)$ time preprocessing, finds an admissible
excess-to-deficit path in $O(r\sqrt{n}\mathop{\mathrm{polylog}}n)$ amortized
time. Summing this cost over all augmentations, we obtain the desired running
time.

###  4.2 An efficient implementation

Recall in the previous sections that we could bound the running time of the
Hungarian search by the size of $\operatorname{supp}(f)$. Here, the number of
active imbalanced nodes at any scale is $O(r)$, and the length of an
augmenting path is also $O(r)$. Therefore one might hope to find an augmenting
path in $O(r\mathop{\mathrm{polylog}}n)$ time, by adapting the algorithms
described in Sections 2 and 3. The challenge is that $\operatorname{supp}(f)$
may have $\Omega(n)$ size, therefore an algorithm which runs in time
proportional to the support size is no longer sufficient. Still, we manage to
implement Hungarian search in time $O(r\sqrt{n}\mathop{\mathrm{polylog}}n)$,
by exploiting a few properties of $\operatorname{supp}(f)$ as described below.

We note that each arc of $\operatorname{supp}(f)$ is admissible with reduced
cost $0$, so we prioritize the relaxation of support arcs as soon as they
arrive in $X\times(V\setminus X)$, over the non-support arcs. This strategy
ensures the following crucial property.

######  Lemma 4.1.

If the support arcs $\operatorname{supp}(f)$ are relaxed as soon as possible,
$\operatorname{supp}(f)$ is acyclic.

Next, similar to Section 3, we call node $u$ _alive_ if (a) $u$ is an active
imbalanced node or (b) if $u$ is incident to an arc of
$\operatorname{supp}(f)$; $u$ is _dead_ otherwise. Unlike in Section 3, once a
node becomes alive it cannot be dead again. Furthermore, a dead node may
become alive only at the beginning of a scale (after the value of $\Delta$ is
reduced). Also, an augmenting path cannot pass through a dead node. Therefore,
we can ignore all dead nodes during Hungarian search, and update the set of
alive/dead nodes at the beginning of a scale.

Let
${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color<EMAIL_ADDRESS>84}}$}}}\subseteq{B}^{\text{\char 96}}$ be the set of nodes that are either
(a) active imbalanced nodes or (b) incident to _exactly one_ arc of
$\operatorname{supp}(f)$. Lemma 4.1 implies that ${B}^{\text{\char
96}}\setminus B^{\text{\char 84}}$ has size $O(r)$. We can therefore find the
min-reduced-cost arc between $X\cap{A}^{\text{\char 96}}$ and
${B}^{\text{\char 96}}\setminus(B^{\text{\char 84}}\cup X)$ using a BCP data
structure as in Section 2, along with lazy potential updates and the rewinding
mechanism. The total time spent by Hungarian search on the nodes of
${B}^{\text{\char 96}}\setminus B^{\text{\char 84}}$ will be
$O(r\mathop{\mathrm{polylog}}n)$. We subsequently focus on handling
$B^{\text{\char 84}}$.

#### Handling $B^{\text{\char 84}}$.

We now describe how we query a min-reduced-cost arc between
$X\cap{A}^{\text{\char 96}}$ and $B^{\text{\char 84}}\setminus X$. Each node
$b\in B^{\text{\char 84}}$ is incident to exactly one arc in
$\operatorname{supp}(f)$. We partition these nodes into clusters depending on
their unique neighbor in $N_{f}$. That is, for a node $a\in{A}^{\text{\char
96}}$, let
${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color<EMAIL_ADDRESS>84}}$}}}\coloneqq\\{b\in B^{\text{\char
84}}\mid{a\mathord{\shortrightarrow}b}\in\operatorname{supp}(f)\\}$. We refer
to $B_{a}^{\text{\char 84}}$ as the _star_ of $a$.

The crucial observation is that $a$ is the only node in $N_{f}$ reachable from
each $b\in B_{a}^{\text{\char 84}}$, so once the Hungarian search reaches a
node of $B_{a}^{\text{\char 84}}$ and thus $a$ (recall we prioritize relaxing
support arcs), the Hungarian search need not visit any other nodes of
$B_{a}^{\text{\char 84}}$, as they will only lead to $a$. Hence, as soon as
one node of $B_{a}^{\text{\char 84}}$ is reached, all other nodes of
$B_{a}^{\text{\char 84}}$ can be discarded from further consideration. Using
this observation, we handle $B^{\text{\char 84}}$ as follows.

We classify each $a\in{A}^{\text{\char 96}}$ as _light_ or _heavy_ : heavy if
$\mathopen{|}B_{a}^{\text{\char 84}}\mathclose{|}\geq\sqrt{n}$, and light if
$\mathopen{|}B_{a}^{\text{\char 84}}\mathclose{|}\leq 2\sqrt{n}$. Note that if
$\sqrt{n}\leq\mathopen{|}B_{a}^{\text{\char 84}}\mathclose{|}\leq 2\sqrt{n}$
then $a$ may be classified as light or heavy. We allow this flexibility to
implement reclassification in a lazy manner. Namely, a light node is
reclassified as heavy once $\mathopen{|}B_{a}^{\text{\char
84}}\mathclose{|}>2\sqrt{n}$, and a heavy node is reclassified as light once
$\mathopen{|}B_{a}^{\text{\char 84}}\mathclose{|}<\sqrt{n}$. This scheme
ensures that the star of $a$ has gone through at least $\sqrt{n}$ updates
between two successive reclassifications, and these updates will pay for the
time spent in updating the data structure when $a$ is re-classified.

For each heavy node $a\in{A}^{\text{\char 96}}\setminus X$, we maintain a BCP
data structure between $B_{a}^{\text{\char 84}}$ and $X\cap{A}^{\text{\char
96}}$. Next, for all light nodes in ${A}^{\text{\char 96}}\setminus X$, we
collect their stars into a single set
${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color<EMAIL_ADDRESS>84}}$}}}\coloneqq\bigcup_{\text{$a$ light}}B_{a}^{\text{\char 84}}$. We
maintain one single BCP data structure between $B_{<}^{\text{\char 84}}$ and
${A}^{\text{\char 96}}\cap X$. Thus, at most $r$ different BCP data structures
are maintained for stars.

Using these data structures, we can compute and relax a min-reduced-cost arc
$v\mathord{\shortrightarrow}w$ between ${A}^{\text{\char 96}}\cap X$ and
$B^{\text{\char 84}}\setminus X$. If $w$ lies in some star $B_{a}^{\text{\char
84}}$, then we also add $a$ into $X$. If $a$ is light, then we delete
$B_{a}^{\text{\char 84}}$ from $B_{<}^{\text{\char 84}}$ and update the BCP
data structure of $B_{<}^{\text{\char 84}}$. If $a$ is heavy, then we stop
querying the BCP data structure of $B_{a}^{\text{\char 84}}$ for the remainder
of the search. Finally, since $a$ becomes part of $X$, $a$ is added to all
$O(r)$ BCP data structures. Recall that $r\leq\sqrt{n}$ by assumption. Adding
arc $v\mathord{\shortrightarrow}w$ thus involves performing $O(\sqrt{n})$
insertion/deletion operations in various BCP data structures, thereby taking
$O(\sqrt{n}\mathop{\mathrm{polylog}}n)$ time.

#### Putting it together.

While proof is omitted, the following lemma bounds the running time of the
Hungarian search.

######  Lemma 4.2.

Assuming all BCP data structures are initialized correctly, the Hungarian
search terminates within $O(r)$ steps, and takes
$O(r\sqrt{n}\mathop{\mathrm{polylog}}n)$ time.

Once an augmenting path is found and the augmentation is performed, the set of
imbalanced nodes and the support arcs change. We thus need to update the sets
$B^{\text{\char 84}}$, $B_{a}^{\text{\char 84}}$s, and $B_{<}^{\text{\char
84}}$. This can be accomplished in $O(r\mathop{\mathrm{polylog}}n)$ amortized
time. When we begin a new Hungarian search, we use the rewinding mechanism to
set various BCP data structures in the right initial state. Finally, when we
move from one scale to another, we also update the sets ${A}^{\text{\char
96}}$ and ${B}^{\text{\char 96}}$. Omitting all the details, we conclude the
following.

######  Lemma 4.3.

Each Hungarian search can be performed in
$O(r\sqrt{n}\mathop{\mathrm{polylog}}n)$ time.

Since there are $O(n\log n)$ augmentations and the flow in the original
network can be recovered from that in the contracted network in
$O(n\mathop{\mathrm{polylog}}n)$ time [3], the total running time of the
algorithm is $O(rn^{3/2}\mathop{\mathrm{polylog}}n)$, as claimed in Theorem
1.3.

#### Acknowledgment.

We thank Haim Kaplan for discussion and suggestion to use Goldberg et al. [9]
algorithm. We thank Debmalya Panigrahi for helpful discussions.

## References

  * [1] Pankaj K. Agarwal, Alon Efrat, and Micha Sharir.  Vertical decomposition of shallow levels in 3-dimensional arrangements and its applications.  _SIAM J. Comput._ 29(3):912–953, 1999.  $\langle$https://doi.org/10.1137/S0097539795295936$\rangle$.
  * [2] Pankaj K. Agarwal, Kyle Fox, Debmalya Panigrahi, Kasturi R. Varadarajan, and Allen Xiao.  Faster algorithms for the geometric transportation problem.  _Proc. 33rd Int. Sympos. Comput. Geom. (SoCG)_ , 7:1–7:16, 2017.  $\langle$https://doi.org/10.4230/LIPIcs.SoCG.2017.7$\rangle$.
  * [3] Pankaj K. Agarwal, Kyle Fox, Debmalya Panigrahi, Kasturi R. Varadarajan, and Allen Xiao.  Faster algorithms for the geometric transportation problem.  Preprint, 2019.  $\langle$http://arxiv.org/abs/1903.08263$\rangle$.
  * [4] Pankaj K. Agarwal and Kasturi R. Varadarajan.  A near-linear constant-factor approximation for Euclidean bipartite matching?  _Proc. 20th Annu. Sympos. Comput. Geom. (SoCG)_ , 247–252, 2004.  $\langle$https://doi.org/10.1145/997817.997856$\rangle$.
  * [5] D. Bertsekas and D. El Baz.  Distributed asynchronous relaxation methods for convex network flow problems.  _SIAM J. Control and Opt._ 25(1):74–85, 1987.  $\langle$https://doi.org/10.1137/0325006$\rangle$.
  * [6] Ran Duan, Seth Pettie, and Hsin-Hao Su.  Scaling algorithms for weighted matching in general graphs.  _ACM Trans. Algorithms_ 14(1):8:1–8:35, 2018.  $\langle$https://doi.org/10.1145/3155301$\rangle$.
  * [7] Shimon Even and Robert E. Tarjan.  Network flow and testing graph connectivity.  _SIAM J. Comput._ 4(4):507–518, 1975.  $\langle$https://doi.org/10.1137/0204043$\rangle$.
  * [8] Harold N. Gabow and Robert E. Tarjan.  Faster scaling algorithms for network problems.  _SIAM J. Comput._ 18(5):1013–1036, 1989.  $\langle$https://doi.org/10.1137/0218069$\rangle$.
  * [9] Andrew V. Goldberg, Sagi Hed, Haim Kaplan, and Robert E. Tarjan.  Minimum-cost flows in unit-capacity networks.  _Theoret. Comput. Sci._ 61(4):987–1010, 2017.  $\langle$https://doi.org/10.1007/s00224-017-9776-7$\rangle$.
  * [10] John E. Hopcroft and Richard M. Karp.  An $n^{5/2}$ algorithm for maximum matchings in bipartite graphs.  _SIAM J. Comput._ 2(4):225–231, 1973.  $\langle$https://doi.org/10.1137/0202019$\rangle$.
  * [11] Haim Kaplan, Wolfgang Mulzer, Liam Roditty, Paul Seiferth, and Micha Sharir.  Dynamic planar Voronoi diagrams for general distance functions and their algorithmic applications.  _Proc. 28th Annu. ACM-SIAM Sympos. Discrete Algorithms (SODA)_ , 2495–2504, 2017.  $\langle$https://doi.org/10.1137/1.9781611974782.165$\rangle$.
  * [12] Andrey Boris Khesin, Aleksandar Nikolov, and Dmitry Paramonov.  Preconditioning for the geometric transportation problem.  Preprint, 2019.  $\langle$http://arxiv.org/abs/1902.08384$\rangle$.
  * [13] Harold W. Kuhn.  The Hungarian method for the assignment problem.  _Naval Research Logistics (NRL)_ 2(1-2):83–97. John Wiley & Sons, 1955.
  * [14] Yin Tat Lee and Aaron Sidford.  Path finding methods for linear programming: Solving linear programs in ${\tilde{O}}(\sqrt{\text{rank}})$ iterations and faster algorithms for maximum flow.  _55th Annu. IEEE Sympos. Found. Comput. Sci. (FOCS)_ , 424–433, 2014\.  $\langle$https://doi.org/10.1109/FOCS.2014.52$\rangle$.
  * [15] Aleksander Mądry.  Navigating central path with electrical flows: From flows to matchings, and back.  _54th Annu. IEEE Sympos. Found. Comput. Sci. (FOCS)_ , 253–262, 2013\.  $\langle$https://doi.org/10.1109/FOCS.2013.35$\rangle$.
  * [16] James B. Orlin.  A faster strongly polynomial minimum cost flow algorithm.  _Operations Research_ 41(2):338–350, 1993.  $\langle$https://doi.org/10.1287/opre.41.2.338$\rangle$.
  * [17] Lyle Ramshaw and Robert E. Tarjan.  A weight-scaling algorithm for min-cost imperfect matchings in bipartite graphs.  _Proc. 53rd Annu. IEEE Sympos. Found. Comput. Sci. (FOCS)_ , 581–590, 2012.  $\langle$https://doi.org/10.1109/FOCS.2012.9$\rangle$.
  * [18] R. Sharathkumar and Pankaj K. Agarwal.  Algorithms for the transportation problem in geometric settings.  _Proc. 23rd Annu. ACM-SIAM Sympos. Discrete Algorithms (SODA)_ , 306–317, 2012.  $\langle$https://dl.acm.org/citation.cfm?id=2095116.2095145$\rangle$.
  * [19] R. Sharathkumar and Pankaj K. Agarwal.  A near-linear time $\epsilon$-approximation algorithm for geometric bipartite matching.  _Proc. 44th Annu. ACM Sympos. Theory Comput. (STOC)_ , 385–394, 2012\.  $\langle$https://doi.org/10.1145/2213977.2214014$\rangle$.
  * [20] Éva Tardos.  A strongly polynomial minimum cost circulation algorithm.  _Combinatorica_ 5(3):247–256, 1985.  $\langle$https://doi.org/10.1007/BF02579369$\rangle$.
  * [21] Pravin M. Vaidya.  Geometry helps in matching.  _SIAM J. Comput._ 18(6):1201–1225, 1989.  $\langle$https://doi.org/10.1137/0218080$\rangle$.
  * [22] Kasturi R. Varadarajan.  A divide-and-conquer algorithm for min-cost perfect matching in the plane.  _Proc. 39th Annu. IEEE Sympos. Found. Comput. Sci. (FOCS)_ , 320–331, 1998.  $\langle$https://doi.org/10.1109/SFCS.1998.743466$\rangle$.


# Spatial Dimensionality Dependence of Heterogeneity, Breakdown of the Stokes-
Einstein Relation and Fragility of a Model Glass-Forming Liquid

Monoj Adhikari  Jawaharlal Nehru Centre for Advanced Scientific Research,
Jakkuru Campus, 560064 Bengaluru, India  Smarajit Karmakar  TIFR Center for
Interdisciplinary Science, Tata Institute of Fundamental Research,36/P
Gopanpally Village, Serilingampally Mandal, RR District, Hyderabad, 500075,
Telangana, India,  Srikanth Sastry<EMAIL_ADDRESS>Jawaharlal Nehru
Centre for Advanced Scientific Research, Jakkuru Campus, 560064, Bengaluru,
India

###### Abstract

We investigate the heterogeneity of dynamics, the breakdown of the Stokes-
Einstein relation and fragility in a model glass forming liquid, a binary
mixture of soft spheres with a harmonic interaction potential, for spatial
dimensions from $3$ to $8$. Dynamical heterogeneity is quantified through the
dynamical susceptibility $\chi_{4}$, and the non-Gaussian parameter
$\alpha_{2}$. We find that the fragility, the degree of breakdown of the
Stokes-Einstein relation, as well as heterogeneity of dynamics, decrease with
increasing spatial dimensionality. We briefly describe the dependence of
fragility on density, and use it to resolve an apparent inconsistency with
previous results.

##  1 Introduction

When a glass-forming liquid is cooled to low temperatures, its dynamics
becomes heterogeneous; spatially correlated clusters of particles move faster
or slower than the average. The growth of dynamic heterogeneity (DH) and its
associated dynamic length scales with the lowering of temperature have been
considered a hallmark of relaxation dynamics in glass forming liquids 1, 2, 3,
4, 5, 6, 7, 8. The role of such heterogeneities in the complex relaxation
dynamics in glass forming liquids, and growing length scales governing these
relaxation processes have been widely investigated 9, 10. An extensively
studied phenomenon, which has been analysed in the context of dynamical
heterogeneity is the breakdown of the Stokes-Einstein relation (SER)
${D\eta\over T}=constant$ (or equivalently, ${D\tau_{\alpha}}=constant$) where
$D$ is the diffusion coefficient, $\eta$ is the shear viscosity, $T$ is the
temperature, and $\tau_{\alpha}$ is the structural relaxation time 11, 12, 13,
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29. Although the
origin of such a breakdown in hopping dynamics have also been investigated 23,
28, 29, significant evidence links the breakdown of the Stokes-Einstein
relation with length scales over which dynamics is heterogeneous23, 30, 29. It
has also been suggested that fragility, which quantifies the degree of non-
Arrhenius increase of relaxation times upon lowering temperature, is also
related to heterogeneous dynamics. Böhmer _et al._ 31 investigated the
correlation between the fragility and heterogeneity of dynamics, by compiling
data for a large number of glass formers. Fragility was quantified by the
fragility index $m$ which measures the steepness of rise of relaxation times
at the glass transition, in an Angell plot32, wherein the logarithm of the
relaxation time is plotted against inverse temperature scaled to the glass
transition temperature ($T_{g}/T$). The KWW exponent, $\beta$, which
characterises stretched exponential relaxation of density fluctuations, was
considered as a measure of the heterogeneity of dynamics. Large values of
fragility index were found to correspond to small values of $\beta$. Such a
correlation between heterogeneity and fragility have been probed in several
works (and also contested 33), typically through consideration of the
relationship between configurational entropy, fragility and cooperative length
scales of dynamics 34, 35, 36, 37. Some of the issues involved have been
addressed within the framework of the random first order transition (RFOT)
theory 38, 39, 40, extensions of mode coupling theory 41, recent exact results
in the limit of infinite spatial dimensions42 and corresponding investigations
of dynamics in variable dimensions 43. In particular, analyses of static and
dynamic behaviour that may be expected in finite dimensions have led to
identification of an upper critical dimension of $d_{u}=8$ above which mean-
field theories provide the correct description44, 45. Within the framework of
the generalised entropy theory as well46, $d=8$ arises as a special dimension,
above which an entropy vanishing transition does not exist at finite
temperature. In this context, it is of interest to understand how aspects of
heterogeneous dynamics, the breakdown of the SER, and fragility depend on
spatial dimensionality. Indeed, some studies have addressed such dependence
47, 25, 27. In particular, Charbonneau _et al._ 27 considered hard sphere
fluids up to $10$ dimensions and showed evidence that the exponent $\omega$ in
the relation $D\sim\tau_{\alpha}^{-1+\omega}$ that quantifies the break down
of the SER vanishes within numerical uncertainty above spatial dimension
$d=8$. However, this remains the only study that has explored the dimension
dependence above $d=4$. A similar study, for a model system with an
interaction potential other than hard core interaction, which permits the
study of both temperature and density dependent behaviour, is therefore
desirable. We undertake such a study in the present work.
We investigate a model glass forming liquid consisting of a binary mixture of
spheres interacting with a harmonic potential, in $3-8$ spatial dimensions. In
the zero-temperature limit, this model has the limiting behaviour of the hard
sphere model whose behaviour is controlled by density alone, while it exhibits
behaviour of dense glass formers at high densities, at finite temperature. We
perform computer simulations and investigate various measures of dynamical
heterogeneity (DH) such as the non-Gaussian parameter, $\alpha_{2}(t)$, the
dynamical susceptibility, $\chi_{4}$ as a function of time for a wide range of
temperatures. We compute the fragilities from the temperature dependence of
the relaxation times, and further investigate the breakdown of the SER from a
comparison of diffusion coefficients and relaxation times. We find a
consistent variation of behaviour as the spatial dimension increases, wherein
the fragility, extent of heterogeneity, and the degree of breakdown of the SER
decrease with increasing spatial dimensionality, consistent with the approach
to mean-field behaviour at $d=8$. We briefly discuss the dependence on density
of fragility and resolve an apparent inconsistency with previously published
results which suggested an increase of fragility with increasing dimension
while the degree of heterogeneity decreased.

The rest of the paper is organized as follows: In section II, we describe the
model and methods related to this study. In section III, we present our main
results, and in section IV, we briefly discuss the density dependence of
fragility and make a comparison with earlier work. Finally, we present a
summary of results and conclusions in section V.

##  2 Simulations details

We investigate a $50:50$ binary mixture of particles that interact with a
harmonic potential given by 48, 49:

$\displaystyle V_{\alpha\beta}(r)$ $\displaystyle=$
$\displaystyle\frac{\epsilon_{\alpha\beta}}{2}\left(1-\frac{r}{\sigma_{\alpha\beta}}\right)^{2},\hskip
44.9554ptr_{\alpha\beta}\leq\sigma_{\alpha\beta}$ (1) $\displaystyle=$
$\displaystyle 0,\hskip 112.10393ptr_{\alpha\beta}>\sigma_{\alpha\beta}$

where $\alpha,\;\beta$ $\in$ (A,B), indicates the type of particle. The two
types of particle differ in their size, with $\sigma_{BB}=1.4\sigma_{AA}$ (and
the diameters are additive), with the interaction strengths being the same for
all pairs. We present results for $3d$ \- $8d$ fixing the density at
$1.3\phi_{J}$, where $\phi_{J}$ is the jamming density. We have used
$\phi_{J}$=$0.645(3$d), $0.467(4$d), $0.319(5$d), $0.209(6$d), $0.133(7$d),
$0.0821(8$d), using estimates by Charbonneau _et al._ 50. In an accompanying
study, we report our estimates of $\phi_{J}$, as well as a dynamical cross
over density $\phi_{0}$. Our estimates of $\phi_{J}$ are close to the values
used here to a high degree of accuracy. Thus, the volume fraction we use in
our simulations are as follows: $0.8384$ ($3$d), $0.6071$ ($4$d), $0.4147$
($5$d), $0.2717$ ($6$d),$0.1729$ ($7$d), $0.1067$ $(8$d). The number density,
$\rho$ is related to the volume fraction $\phi$ for the binary mixture in the
following way

$\phi=\rho
2^{-d}\frac{\pi^{d/2}}{\Gamma(1+\frac{d}{2})}((c_{A}\sigma^{d}_{AA}+c_{B}\sigma^{d}_{BB})$
(2)

where $\rho=N/V$, with $N$ being the number of particles, and $V$ the volume,
and the fractions $c_{A}=c_{B}=1/2$. The corresponding number densities are
following: $0.8556$ ($3$d),$0.8132$ ($4$d), $0.7904$ ($5$d), $0.7891$ ($6$d),
$0.8114$ ($7$d), $0.8554$ ($8$d).

The system size is fixed at $5000$ particles, which is large enough that the
linear dimension $L$ is $>2\sigma_{BB}$ in all dimensions. Molecular dynamics
(MD) simulations are performed in a cubic box with periodic boundary
conditions in the constant number, volume, and temperature (NVT) ensemble. The
integration time step was fixed at $dt=0.01$. Temperatures are kept constant
using the Brown and Clarke 51 algorithm. The data, presented here, have run
lengths of around $100\tau$ (where $\tau$ is the relaxation time, defined
below). We present results that are averaged over five independent samples.
For results over a range of densities which we discuss in section IV, results
are from 1-2 independent samples at densities other than those mentioned
above. We use reduced units with the small particle diameter, $\sigma_{AA}$,
as the unit of length, $\epsilon_{AA}$ as the energy unit, and
$\sqrt{\frac{\sigma_{AA}^{2}m_{AA}}{\epsilon_{AA}}}$ as time unit, where
$m_{AA}$ is the mass which is set to unity.

##  3 Results

###  3.1 Fragility in different dimensions

We quantify the microscopic dynamics by computing the overlap function, which
is defined (for the $B$ particles) by:

$q(t)=\frac{1}{N_{B}}\sum_{i}^{N_{B}}w(|\textbf{r}_{i}(t_{0})-\textbf{r}_{i}(t+t_{0})|)\quad\mbox{where}\quad
w(x)=\begin{cases}1.0&\text{if $x\leq a$}\\\ 0&\text{otherwise.}\end{cases}\\\
$

Here a is the cut-off within which particle positions are treated as
indistinguishable52. We choose the parameter a in such a way that $a^{2}$ is
close to the plateau value of the mean squared displacement (MSD). The choice
of a is further refined by considering the behaviour of $\chi_{4}$ (defined
below), to identify the value of a for which the peak value of $\chi_{4}$ is
maximum with the choice of a. Fig. 1 illustrates the choice for three
dimensions (3d), and the a value for other dimensions is chosen by a similar
procedure. We choose the parameter values $a=0.48,0.50,0.55,0.60,0.75$ ,and
$0.80$ for $3$d, $4$d, $5$d, $6$d, $7$d, and $8$d, respectively.

Figure 1: Choice of the cut-off parameter a for calculating the overlap
function, illustrated for $3$d, with $T=0.0038$ as the representative
temperature. Left: We show the mean squared displacement (MSD) as a function
of time. The horizontal red line corresponds to $a^{2}=0.2304$. Right:
$\chi_{4}(t)$ is shown as a function of time for different a values. The
choice $a=0.48$ leads to the maximum peak value of $\chi_{4}$.

Figure 2: Relaxation time as a function of temperature for various dimensions.
Left: Points (dot) represent simulation data along with fits (lines) to the
VFT form. Right: Angell plot, of logarithm of the relaxation times, plotted
against $T/T_{g}$ where the glass transition temperature $T_{g}$ is chosen to
be the temperature at which the relaxation time reaches the value $10^{6}$.
The fragility is highest for $3d$, with the fragility decreasing with
increasing spatial dimensionality.

Figure 3: Left: Kinetic fragility is plotted as a function of spatial
dimensionality at fixed density $\phi=1.3\phi_{J}$, which decreases with
increasing dimensionality. Right: The divergence temperature $T_{VFT}$ is
plotted as a function of spatial dimensionality at fixed density
$\phi=1.3\phi_{J}$ which also decreases with increasing dimensionality.

We calculate the relaxation times by considering the overlap function for the
$B$ particles. The relaxation time, $\tau_{\alpha}$ is computed as the time at
which $\langle q(t)\rangle=1/e$, where $\langle\cdots\rangle$ refers to an
ensemble average (we average over initial times and over samples). We compute
the relaxation times for a wide range of temperatures in each dimension. The
relaxation times exhibit super-Arrhenius temperature dependence, the strength
of which is quantified by the kinetic fragility. Here, we estimate the kinetic
fragility from Vogel-Fulcher-Tammann (VFT) fits to the temperature dependence
of the $\alpha$ relaxation times:

$\tau_{\alpha}=\tau_{0}\exp\left[\frac{1}{K_{VFT}\left(\frac{T}{T_{VFT}}-1\right)}\right],$
(3)

where $K_{VFT}$ is the kinetic fragility of the system and $T_{VFT}$ is the
temperature at which the relaxation time diverges by extrapolation. In Fig. 2,
we show relaxation time as a function of temperature in a semi-log plot for
each dimension, against temperature, as well as against scaled inverse
temperature $T_{g}/T$, in an Angell plot. The glass transition temperature
$T_{g}$ is defined as the temperature where relaxation time becomes $10^{6}$.
The pre-factor in Eq. 3, $\tau_{0}\approx 10$ in all dimensions, and thus its
inclusion or otherwise in defining $T_{g}$ does not alter the observed
behavior. From the Angell plot, it is apparent that the liquid becomes more
fragile as the spatial dimension increases. The kinetic fragility $K_{VFT}$
and the divergence temperature $T_{VFT}$ are plotted as a function of
dimension in Fig. 3, which shows that both $K_{VFT}$ and $T_{VFT}$ are
decreasing functions of spatial dimensionality.

###  3.2 Heterogeneity in dynamics

We next investigate the heterogeneity in dynamics by two different measures of
heterogeneity: 1. The dynamical susceptibility, $\chi_{4}$, and 2. The non-
Gaussian parameter, $\alpha_{2}$.

Figure 4: The dynamical susceptibility $\chi_{4}(t)$ is plotted as a function
of time for different temperatures for spatial dimensions $3-8$.

####  3.2.1 Dynamical susceptibility, $\chi_{4}$

Dynamical susceptibility, $\chi_{4}$, which measures the fluctuations in the
overlap function $q(t)$, is defined by:

$\chi_{4}(t)=N_{B}\left[\langle q(t)^{2}\rangle-\langle
q(t)\rangle^{2}\right]$ (4)

where the average is over initial configurations and the independent samples.

Figure 5: Left: The peak value of $\chi_{4}(t)$, $\chi_{4}^{peak}$, is shown
as a function of relaxation time, $\tau_{\alpha}$. Inset: By scaling
$\chi_{4}^{peak}$ values, the data for different dimensions are collapsed onto
a master curve. A power law fit (red line) provides a reasonable description
for most of the temperature ($\tau_{\alpha}$) range, with exponent $0.47$.
Right: The time at which $\chi_{4}(t)$ is maximum, $t_{\chi_{4}}^{max}$ is
plotted as a function of $\tau_{\alpha}$ for different dimensions. The data
for different dimensions overlap, and demonstrate that
$t_{\chi_{4}}^{max}\sim\tau_{\alpha}$.

As it has been demonstrated extensively, the time dependence of $\chi_{4}(t)$
is non-monotonic, and exhibits a peak value, ($\chi_{4}^{peak}$), at a time
that is proportional to the alpha relaxation time. In Fig. 4, we show
$\chi_{4}(t)$ against time for different temperatures in each dimension. The
peak value of $\chi_{4}^{peak}$ as well as the time at which it occurs,
$t_{\chi_{4}}^{max}$, increase strongly upon a decrease in temperature,
indicating that the heterogeneity of dynamics increases with a decrease in
temperature, and is maximum at a time scale that increases in proportion to
$\tau_{\alpha}$. To compare the degree of heterogeneity for different spatial
dimensions, we show, in Fig. 5 (Left panel), $\chi_{4}^{peak}$ as a function
of $\tau_{\alpha}$ for each dimension. For a given $\tau_{\alpha}$,
$\chi_{4}^{peak}$ decreases with increasing spatial dimension implying that
heterogeneity decreases with increasing spatial dimensionality. We also
observe that $\chi_{4}^{peak}$ shows a power law dependence on $\tau_{\alpha}$
at higher temperatures as $\chi_{4}^{peak}\sim\tau_{\alpha}^{z}$, with $z$
being the power-law exponent. Deviations from power law behaviour is mostly
observed at lower temperatures. This behaviour is consistent with previous
observations in three dimensions 53, 54, 8, 52, including polymeric glass
formers 55 and exponent $z$ is found to be close to $0.62$ for binary hard
sphere fluids 54, 8 whereas $z=0.51$ has recently been found for a model
liquid that aims to tune the degree of mean field character 56. Our estimate
is on the lower side of these values, but close to that reported in 56. Since
the power law regime is limited in extent, there is room for error in the
exact determination of the exponent. Remarkably, however, we find that the
exponent of the power law, $z$, is the same in all dimensions for our studied
model, as evidenced by the data collapse, regardless of the precise value. The
expectation of a power law dependence arises, for example, from inhomogeneous
mode coupling theory 41, and the deviations from the power law are understood
to be a consequence of the role played by activated processes at low
temperatures. Thus, the observation of a common exponent describing the power
law dependence of $\chi_{4}^{peak}$ on $\tau_{\alpha}$ should perhaps be seen
as mean-field behaviour that does not depend on spatial dimensions and thus
not surprising. Nevertheless, to our knowledge, such a universal behaviour has
not previously been reported across the range of spatial dimensions that we
investigate.

####  3.2.2 Non-Gaussian parameter, $\alpha_{2}(t)$

Figure 6: The non-Gaussian parameter $\alpha_{2}$ is plotted against time for
different temperatures and spatial dimensions $3-8$.

Figure 7: Left: The peak value of $\alpha_{2}(t)$, $\alpha_{2}^{peak}$, is
shown as a function of relaxation time, $\tau_{\alpha}$ for different
dimensions. Inset: By scaling $\alpha_{2}^{peak}$ values, the data for
different dimensions are collapsed onto a master curve. A power law fit (red
line) provides a reasonable description for most of the temperature
($\tau_{\alpha}$) range, with exponent $0.43$. Right: The time at which
$\alpha_{2}(t)$ is maximum, $t^{*}$ is plotted as a function of $(D/T)^{-1}$
for different dimensions. The data for different dimensions overlap, and
demonstrate that $t^{*}\sim(D/T)^{-1}$, albeit with small deviations apparent
at low temperatures.

Next, we investigate another measure of heterogeneity, the non-Gaussian
parameter, $\alpha_{2}(t)$, for different dimensions. As previously discussed
in detail 36, 57, 55, $\alpha_{2}(t)$ and $\chi_{4}$ correspond to distinct
aspects of heterogeneity, associated with correlated clusters of mobile, and
immobile, particles respectively.. The non-Gaussian parameter, $\alpha_{2}(t)$
measures the deviation of the van Hove distribution of displacements of the
particle in time $t$ from Gaussian form, expected for spatially homogeneous
dynamics, and is given by

$\displaystyle\alpha_{2}(t)$ $\displaystyle=$ $\displaystyle
C_{d}\frac{\langle r^{4}(t)\rangle}{\langle r^{2}(t)\rangle^{2}}-1$ (5)
$\displaystyle\langle r^{2n}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{N}\langle(\vec{r}_{i}(t)-\vec{r}_{i}(0))^{2n}\rangle$
(6)

where $C_{d}$ is a spatial dimension dependent coefficient to ensure that
$\alpha_{2}(t)=0$ when the distribution of displacements is a Gaussian.
Similar to $\chi_{4}$, $\alpha_{2}$ also shows non-monotonic behaviour with
respect to time. However, the characteristic time $t^{*}$ at which
$\alpha_{2}(t)$ is maximum is smaller than $\tau_{\alpha}$, and has been
demonstrated to be proportional to a time scale determined by the diffusion
coefficient, $(D/T)^{-1}$ 36. In Fig. 6, we show $\alpha_{2}(t)$ against time
for different temperatures and spatial dimensions $3-8$. We see that
$\alpha_{2}^{peak}$ increases with a decrease in temperature for all spatial
dimensions. Similarly to $\chi_{4}$, we report $\alpha_{2}^{peak}$ against
$\tau_{\alpha}$ for spatial dimension $3-8$ in Fig. 7 (Left panel). We see
that for a given $\tau_{\alpha}$, $\alpha_{2}^{peak}$ also decreases with
increasing dimensionality. This again implies that heterogeneity decreases
with increasing spatial dimensionality. Similar to $\chi_{4}^{peak}$,
$\alpha_{2}^{peak}$ also displays a power law dependence on $\tau_{\alpha}$ at
higher dimensions, with deviations at lower temperatures. The inset of Fig. 7
(Left panel) demonstrates that for $\alpha_{2}^{peak}$, the exponent of the
power law is $0.43$ and is a good description of the data for all dimensions.
In 58 the behaviour of $\alpha_{2}^{peak}$ vs. $\tau_{\alpha}$ was fitted to
two power laws, with exponents $0.8$ (high temperatures) and $0.3$ (low
temperatures). While we find the exponent to be $0.43$ convincingly over two
decades of (high to moderate temperature) relaxation times, we do find that an
exponent of $0.3$ is a good description of low temperature data. In Fig. 7
(Right panel), we show the time $t^{*}$ against $(D/T)^{-1}$, where $D$ is
diffusivity, for different spatial dimensions, confirming the validity of the
relation $t^{*}\sim(D/T)^{-1}$ beyond three dimensions 36. The observed
relationship between $t^{*}$ and $D/T$ has been found to be valid in many
different glass formers as well as other materials 36, 25, 59, 60, and our
results show that it is valid in different dimensions as well. We note,
however, that an exponent other than $-1$ has been reported recently for a
metallic glass former 61.

###  3.3 The Breakdown of the Stokes-Einstein Relation

A much studied phenomenon associated with glassy behavoiur is the violation or
breakdown of the Stokes-Einstein relation (SER), which relates the
translational diffusion coefficient ($D$) of a Brownian particle to the shear
viscosity $\eta$ of the surrounding liquid at a temperature T: $D=mk_{B}T/c\pi
R\eta$, where $m$ is the mass and $R$ is the radius of the particle, $T$ is
the temperature of the liquid, and the factor $c$ is a constant which depends
on the boundary condition at the surface of the Brownian particle. It is been
observed in several investigations that the SER is also satisfied when one
considers the self-diffusion of particles in a liquid at relatively high
temperatures (The caveats and the extent to which such a statement is valid
have also been discussed, e. g. 27). However, as temperature is decreased
towards the glass transition, the SER is observed to break down. As mentioned
in the introduction, violations of the SER, which can be expressed as
${D\eta\over T}=constant$, have been investigated considering $\tau_{\alpha}$
in place of ${\eta\over T}$, expressing the SER as $D\tau_{\alpha}=constant$
25, 27, 29. This equivalence has been validated by computing the viscosity
$\eta$ and comparing with $\tau_{\alpha}$, at the wave vector corresponding to
the peak of the structure factor 25. Further, several works have considered
relaxation behaviour as a function of the wave vector $k$ 20, 25, 27, 29,
either through the $k$ dependent viscosity 20 or relaxation times computed as
a function of $k$. In 29, it was shown that for a given $k$, violation of the
SER arise when $k^{-1}$ falls below a length scale characterising dynamical
heterogeneity. In the present work, we do not investigate the $k$ dependence
of the violation of SER, but consider only the $\tau_{\alpha}$ defined above,
which empirically is equivalent to considering the SER between diffusion
coefficients and viscosity, as mentioned above. The enhancement of the
diffusion coefficient was obtained by Kim and Keyes 20 through a mode coupling
expression for the diffusion coefficient as an integral over the inverse of
the $k$ dependent viscosity $\eta(k)$. Thus, the reduction of $\eta(k)$ with
respect to the hydrodynamic value is offered as a compelling explanation of
the violation of the SER. Indeed, intuitively, the observations employing
relaxation times obtained as a function of wave-vector 25, 27, 29 are
consistent with such an explanation. An investigation of the relationship
between these different approaches has not, however, been systematically
carried out, and would be interesting to perform. The breakdown of SER is
characterized by an exponent $\omega$ that describes a scaling
$D\sim\tau_{\alpha}^{-1+\omega}$. As mentioned above, a limited number of
previous studies 47, 25, 27 have considered SER and the breakdown thereof as a
function of spatial dimension. Charbonneau _et al._ 27 have performed a
hydrodynamic analysis of SER for varying spatial dimension, as well as
numerical investigations up to $d=10$ for hard sphere liquids. Here, we
examine validity or breakdown of the SER employing the self diffusion
coefficients $D$ and the $\tau_{\alpha}$ described above, as a function of
spatial dimensionality.

Figure 8: $D\tau_{\alpha}$ is plotted against $\tau_{\alpha}$ in a log-log
plot. The low temperature data is fitted with the form:
$D\tau_{\alpha}\sim\tau_{\alpha}^{\omega}$. From the fit, we obtain $\omega$
for each spatial dimension.

In Fig. 8, we show the diffusivity, $D$ multiplied by the relaxation time,
$\tau_{\alpha}$, against $\tau_{\alpha}$ in log-log plot. We observe that
$D\tau_{\alpha}$ is roughly constant for small $\tau_{\alpha}$ (high
temperature), at least for $3d,4d$ and $5d$, but display power law behaviour
with a finite $\omega$ for large $\tau_{\alpha}$. We obtain the exponent
$\omega$ by power-law fits of the form
$D\tau_{\alpha}\sim\tau_{\alpha}^{\omega}$ for the data in the low temperature
regime. In Fig. 9 (Left panel), we show $\omega$ as a function of spatial
dimension. We find that $\omega$ is large for $3$d, and decreases with
increasing spatial dimensionality following a relation, $\omega\sim(d-d_{c})$.
We see $\omega$ becomes zero by extrapolation of the linear form at d$=8$,
consistently with the idea that the upper critical dimension $d_{u}=8$, and
consistently with results for the hard sphere fluid 62. The numerical values
of $\omega$ show a deviation from the linear fit for $d=7,8$ and do not vanish
for $d=8$, as also seen in 62. With the results available, we cannot further
probe this issue. Improved numerical results and performing simulations at
higher dimensions than $8$ will permit more precise statements in this regard.
In Fig. 9 (Right panel), we show the breakdown temperature for SER relation,
$T_{SEB}$, as a function of spatial dimension. $T_{SEB}$ is defined as the
temperature below which $D\tau_{\alpha}$ exceeds the high temperature value by
$7.5\%$..

Figure 9: Left panel: The exponent $\omega$ is plotted against the spatial
dimensions. The exponent $\omega$ decreases with increasing spatial
dimensionality. Right panel: $T_{SEB}$, the temperature where SER breaks down
is plotted against spatial dimension. $T_{SEB}$ also decreases with increasing
spatial dimensionality.

###  3.4 Density Dependence and Comparison with Previous Results

We next briefly consider the dependence of the fragility on the density at
which the liquids are studied, and address an apparent inconsistency with
previous results. Further details of density dependence are discussed in
detail in an accompanying paper. In Fig. 10, we show the kinetic fragility
$K_{VFT}$ and the divergence temperature $T_{VFT}$ against density, scaled
with $\phi_{J}$. We note that the kinetic agility $K_{VFT}$ decreases and
nearly vanishes as the density is decreased towards $\phi_{J}$ (Indeed, for
higher spatial dimensions, such vanishing appears to occur for densities
higher than $\phi_{J}$, whose significance is discussed elsewhere), while at
any fixed density, the fragility is a decreasing function of spatial
dimensionality, consistently with the results discussed already for
$\phi=1.3\phi_{J}$. Similarly to $K_{VFT}$, $T_{VFT}$ also decreases as the
density is lowered, while being smaller for higher dimensions at fixed
density. Thus, in comparing behaviour as a function of dimensionality, care
must be exercised to compare results at the same scaled densities.

Figure 10: Left panel: The kinetic fragility $K_{VFT}$ is shown as a function
of scaled density $\phi/\phi_{J}$ for spatial dimensions $3-8$. Right panel:
$T_{VFT}$ is shown as a function of scaled density $\phi/\phi_{J}$ for spatial
dimensions $3-8$.

The observed dependence on density and spatial dimensionality helps explain an
apparent inconsistency with results discussed by Sengupta _et al._25. In 25,
simulation results were shown for the Kob-Andersen (KA) binary Lennard-Jones
mixture, and it was observed that liquids in $4d$ were less heterogeneous than
in $3d$ (consistently with results here), but had larger fragility than in
$3d$, which is not consistent with the present observations that the fragility
too decreases with increasing spatial dimensionality.

Figure 11: Left panel: Relaxation times are shown in an Angell plot for (a)
$3$d at two different densities, $1.14\phi_{J}$, $1.3\phi_{J}$ and (b) for
$4$d at density $1.3\phi_{J}$. The kinetic fragility at $1.14\phi_{J}$ ($3$d)
is lower than the $3d$ and $4d$ systems at density $1.3\phi_{J}$. Right panel:
$\chi_{4}^{peak}$ is plotted against $\tau_{\alpha}$ for $3$d at two different
densities, $1.14\phi_{J}$ and $1.3\phi_{J}$ and for $4$d at $1.3\phi_{J}$.
$\chi_{4}^{peak}$ values for $4$d at $1.3\phi_{J}$ are lower than for the $3d$
systems at both $1.14\phi_{J}$ and $1.3\phi_{J}$.

As noted above, the fragility as well as the heterogeneity depends upon the
density for a given spatial dimension, and thus, to compare results in
different dimensions, one must consider appropriate densities. For the $KA$
system, we do not have a jamming density to provide an appropriate scale, and
instead, we use the zero temperature limit of the spinodal density for
reference 63. The spinodal density for the KA system is $\sim 1.1$ for $3$d
whereas it is $\sim 1.4$ for $4$d. The simulations in 25 were performed at
higher densities, $1.2$ for $3$d which is $1.09$ times the spinodal density
whereas for $4$d the density of $1.6$ was employed, which is $1.14$ times the
spinodal density. Thus, the scaled density in $4d$ is higher than the scaled
density in $3d$, which leads to a higher fragility in the higher dimension. To
illustrate this possibility, we consider for the soft sphere system a lower
density of $1.14\phi_{J}$ in $3d$, in addition to $1.3\phi_{J}$, and compare
with results in $4d$ at $1.3\phi_{J}$.

In Fig. 11 (Left panel), we show an Angell plot for three different cases:
$1.14\phi_{J}$ and $1.30\phi_{J}$ at $3$d and $1.30\phi_{J}$ at $4$d. We note
that the $3d$ system at $1.14\phi_{J}$ has a lower fragility (as confirmed by
calculating $K_{VFT}$) than the $4d$ system at $1.3\phi_{J}$. Plotting
$\chi_{4}^{peak}$ against $\tau_{\alpha}$ for the same three cases (Right
panel, Fig. 11), we see that the $4d$ system has the lowest heterogeneity.
Thus, comparing the $3d$ system at $1.14\phi_{J}$ with the $4d$ system at
$1.3\phi_{J}$ would lead to the conclusion that the heterogeneity in $3d$ is
higher than in $4d$ while the fragility is higher in $4d$, whereas comparison
at the same scaled density would lead to the conclusion that both the
fragility and heterogeneity would decrease with increasing spatial
dimensionality. We therefore conclude that the results in 25 can be understood
consistently with our present results by noting the choice of densities in 25.

##  4 Summary and Conclusions

In conclusion, we have investigated the relationship between fragility,
heterogeneity and the breakdown of the Stokes-Einstein relation in different
spatial dimensions. Our results show that at fixed density, the fragility, the
degree of heterogeneity and the degree of violation of the Stokes-Einstein
relation decrease with increasing spatial dimensionality. The heterogeneity
measures $\chi_{4}^{peak}$ and $\alpha_{2}^{peak}$ depend on the relaxation
time $\tau_{\alpha}$ at high and moderate temperatures in a power law fashion,
with power law exponents that do not depend on spatial dimensionality. The
exponent $\omega$ that characterises the breakdown of the Stokes-Einstein
relationship displays a nearly linear relationship with spatial dimensions
that corresponds to a vanishing of $\omega$ at $d=8$, consistently with the
idea that $d=8$ represents the upper critical dimension. The $\omega$ values
in $d=7,8$ display small deviations from such linear behaviour with $d$, which
requires further investigation including studies in dimensions above $8$. We
show that fragilities decrease with density at all spatial dimensions,
approaching Arrhenius behaviour close to the jamming density. The observed
density dependence helps rationalise earlier results that suggested that
fragility and heterogeneity may vary in opposite ways as a function of spatial
dimensionality.

We acknowledge the Thematic Unit of Excellence on Computational Materials
Science, and the National Supercomputing Mission facility (Param Yukti) at the
Jawaharlal Nehru Center for Advanced Scientific Research for computational
resources. SK acknowledges support from Swarna Jayanti Fellowship grants
DST/SJF/PSA-01/2018-19 and SB/SFJ/2019-20/05. SS acknowledges support through
the JC Bose Fellowship (JBR/2020/000015) SERB, DST (India).

## References

  * Sillescu 1999 Sillescu, H. _Journal of Non-Crystalline Solids_ 1999, _243_ , 81–108
  * Ediger 2000 Ediger, M. D. _Annual review of physical chemistry_ 2000, _51_ , 99–128
  * Berthier et al. 2011 Berthier, L.; Biroli, G.; Bouchaud, J.-P.; Cipelletti, L.; van Saarloos, W. _Dynamical heterogeneities in glasses, colloids, and granular media_ ; OUP Oxford, 2011; Vol. 150
  * Kob et al. 1997 Kob, W.; Donati, C.; Plimpton, S. J.; Poole, P. H.; Glotzer, S. C. _Physical review letters_ 1997, _79_ , 2827
  * Perera and Harrowell 1998 Perera, D. N.; Harrowell, P. _Physical review letters_ 1998, _81_ , 120
  * Yamamoto and Onuki 1998 Yamamoto, R.; Onuki, A. _Physical review letters_ 1998, _81_ , 4915
  * Donati et al. 2002 Donati, C.; Franz, S.; Glotzer, S. C.; Parisi, G. _Journal of non-crystalline solids_ 2002, _307_ , 215–224
  * Flenner et al. 2014 Flenner, E.; Staley, H.; Szamel, G. _Physical review letters_ 2014, _112_ , 097801
  * Karmakar et al. 2014 Karmakar, S.; Dasgupta, C.; Sastry, S. _Annu. Rev. Condens. Matter Phys._ 2014, _5_ , 255–284
  * Karmakar et al. 2016 Karmakar, S.; Dasgupta, C.; Sastry, S. _Reports on Progress in Physics_ 2016, _79_ , 016601
  * Rössler 1990 Rössler, E. _Phys. Rev. Lett._ 1990, _65_ , 1595–1598
  * Fujara et al. 1992 Fujara, F.; Geil, B.; Sillescu, H.; Fleischer, G. _Zeitschrift für Physik B Condensed Matter_ 1992, _88_ , 195–204
  * Thirumalai and Mountain 1993 Thirumalai, D.; Mountain, R. D. _Physical Review E_ 1993, _47_ , 479
  * Stillinger and Hodgdon 1994 Stillinger, F. H.; Hodgdon, J. A. _Physical review E_ 1994, _50_ , 2064
  * Tarjus and Kivelson 1995 Tarjus, G.; Kivelson, D. _The Journal of chemical physics_ 1995, _103_ , 3071–3073
  * Cicerone and Ediger 1996 Cicerone, M. T.; Ediger, M. D. _The Journal of chemical physics_ 1996, _104_ , 7210–7218
  * Berthier et al. 2004 Berthier, L.; Chandler, D.; Garrahan, J. P. _EPL (Europhysics Letters)_ 2004, _69_ , 320
  * Berthier 2004 Berthier, L. _Physical Review E_ 2004, _69_ , 020201
  * Jung et al. 2004 Jung, Y.; Garrahan, J. P.; Chandler, D. _Phys. Rev. E_ 2004, _69_ , 061205
  * Kim and Keyes 2005 Kim, J.; Keyes, T. _The Journal of Physical Chemistry B_ 2005, _109_ , 21445–21448
  * Kumar et al. 2006 Kumar, S. K.; Szamel, G.; Douglas, J. F. _The Journal of chemical physics_ 2006, _124_ , 214501
  * Becker et al. 2006 Becker, S. R.; Poole, P. H.; Starr, F. W. _Physical review letters_ 2006, _97_ , 055901
  * Chong 2008 Chong, S.-H. _Physical Review E_ 2008, _78_ , 041501
  * Chong and Kob 2009 Chong, S.-H.; Kob, W. _Physical review letters_ 2009, _102_ , 025702
  * Sengupta et al. 2013 Sengupta, S.; Karmakar, S.; Dasgupta, C.; Sastry, S. _The Journal of chemical physics_ 2013, _138_ , 12A548
  * Sengupta and Karmakar 2014 Sengupta, S.; Karmakar, S. _The Journal of chemical physics_ 2014, _140_ , 224505
  * Charbonneau et al. 2013 Charbonneau, B.; Charbonneau, P.; Jin, Y.; Parisi, G.; Zamponi, F. _The Journal of chemical physics_ 2013, _139_ , 164502
  * Charbonneau et al. 2014 Charbonneau, P.; Jin, Y.; Parisi, G.; Zamponi, F. _Proceedings of the National Academy of Sciences_ 2014, _111_ , 15025–15030
  * Parmar et al. 2017 Parmar, A. D.; Sengupta, S.; Sastry, S. _Physical review letters_ 2017, _119_ , 056001
  * Nandi and Bhattacharyya 2019 Nandi, M. K.; Bhattacharyya, S. M. _Journal of Physics: Condensed Matter_ 2019, _32_ , 064001
  * Böhmer et al. 1993 Böhmer, R.; Ngai, K.; Angell, C. A.; Plazek, D. _The Journal of chemical physics_ 1993, _99_ , 4201–4209
  * Angell 1991 Angell, C. _Journal of Non-Crystalline Solids_ 1991, _131_ , 13–31
  * Nielsen et al. 2009 Nielsen, A. I.; Christensen, T.; Jakobsen, B.; Niss, K.; Olsen, N. B.; Richert, R.; Dyre, J. C. _The Journal of chemical physics_ 2009, _130_ , 154508
  * Douglas et al. 2006 Douglas, J. F.; Dudowicz, J.; Freed, K. F. _The Journal of chemical physics_ 2006, _125_ , 144907
  * Dudowicz et al. 2007 Dudowicz, J.; Freed, K. F.; Douglas, J. F. _Advances in Chemical Physics_ ; John Wiley & Sons, Ltd, 2007; Chapter 3, pp 125–222
  * Starr et al. 2013 Starr, F. W.; Douglas, J. F.; Sastry, S. _The Journal of chemical physics_ 2013, _138_ , 12A541
  * Betancourt et al. 2013 Betancourt, B. A. P.; Douglas, J. F.; Starr, F. W. _Soft Matter_ 2013, _9_ , 241–254
  * Kirkpatrick and Thirumalai 1987 Kirkpatrick, T. R.; Thirumalai, D. _Physical review letters_ 1987, _58_ , 2091
  * Kirkpatrick and Thirumalai 1988 Kirkpatrick, T. R.; Thirumalai, D. _Physical Review A_ 1988, _37_ , 4439
  * Kirkpatrick et al. 1989 Kirkpatrick, T. R.; Thirumalai, D.; Wolynes, P. G. _Physical Review A_ 1989, _40_ , 1045
  * Biroli et al. 2006 Biroli, G.; Bouchaud, J.-P.; Miyazaki, K.; Reichman, D. R. _Physical review letters_ 2006, _97_ , 195701
  * Charbonneau et al. 2017 Charbonneau, P.; Kurchan, J.; Parisi, G.; Urbani, P.; Zamponi, F. _Annual Review of Condensed Matter Physics_ 2017, _8_ , 265–288
  * Manacorda et al. 2020 Manacorda, A.; Schehr, G.; Zamponi, F. _The Journal of Chemical Physics_ 2020, _152_ , 164506
  * Biroli and Bouchaud 2007 Biroli, G.; Bouchaud, J.-P. _Journal of Physics: Condensed Matter_ 2007, _19_ , 205101
  * Franz et al. 2012 Franz, S.; Jacquin, H.; Parisi, G.; Urbani, P.; Zamponi, F. _Proceedings of the National Academy of Sciences_ 2012, _109_ , 18725–18730
  * Xu et al. 2016 Xu, W.-S.; Douglas, J. F.; Freed, K. F. _Adv. Chem. Phys_ 2016, _161_ , 443–497
  * Eaves and Reichman 2009 Eaves, J. D.; Reichman, D. R. _Proceedings of the National Academy of Sciences_ 2009, _106_ , 15171–15175
  * Durian 1995 Durian, D. J. _Physical review letters_ 1995, _75_ , 4780
  * Berthier and Witten 2009 Berthier, L.; Witten, T. A. _EPL (Europhysics Letters)_ 2009, _86_ , 10001
  * Charbonneau et al. 2011 Charbonneau, P.; Ikeda, A.; Parisi, G.; Zamponi, F. _Physical review letters_ 2011, _107_ , 185702
  * Brown and Clarke 1984 Brown, D.; Clarke, J. _Molecular Physics_ 1984, _51_ , 1243–1252
  * Lačević et al. 2003 Lačević, N.; Starr, F. W.; Schrøder, T.; Glotzer, S. _The Journal of chemical physics_ 2003, _119_ , 7372–7387
  * Karmakar et al. 2009 Karmakar, S.; Dasgupta, C.; Sastry, S. _Proceedings of the National Academy of Sciences_ 2009, _106_ , 3675–3679
  * Flenner and Szamel 2010 Flenner, E.; Szamel, G. _Physical review letters_ 2010, _105_ , 217801
  * Xu et al. 2020 Xu, W.-S.; Douglas, J. F.; Xu, X. _Macromolecules_ 2020, _53_ , 4796–4809
  * Nandi et al. 2021 Nandi, U. K.; Kob, W.; Maitra Bhattacharyya, S. _The Journal of Chemical Physics_ 2021, _154_ , 094506
  * Xu et al. 2016 Xu, W.-S.; Douglas, J. F.; Freed, K. F. _Macromolecules_ 2016, _49_ , 8355–8370
  * Wang et al. 2018 Wang, L.; Xu, N.; Wang, W.; Guan, P. _Physical review letters_ 2018, _120_ , 125502
  * Wang et al. 2019 Wang, X.; Xu, W.-S.; Zhang, H.; Douglas, J. F. _The Journal of chemical physics_ 2019, _151_ , 184503
  * Zhang et al. 2019 Zhang, H.; Wang, X.; Chremos, A.; Douglas, J. F. _The Journal of chemical physics_ 2019, _150_ , 174506
  * Zhang et al. 2021 Zhang, H.; Wang, X.; Yu, H.-B.; Douglas, J. F. _The European Physical Journal E_ 2021, _44_ , 1–30
  * Charbonneau et al. 2012 Charbonneau, P.; Ikeda, A.; Parisi, G.; Zamponi, F. _Proceedings of the National Academy of Sciences_ 2012, _109_ , 13939–13943
  * Sastry 2000 Sastry, S. _Physical Review Letters_ 2000, _85_ , 590