# Signatures of Ultrafast Reversal of Excitonic Order in Ta2NiSe5

H. Ning  Institute for Quantum Information and Matter, California Institute of
Technology, Pasadena, CA 91125  Department of Physics, California Institute of
Technology, Pasadena, CA 91125  O. Mehio  Institute for Quantum Information
and Matter, California Institute of Technology, Pasadena, CA 91125  Department
of Physics, California Institute of Technology, Pasadena, CA 91125  M.
Buchhold  Institute for Quantum Information and Matter, California Institute
of Technology, Pasadena, CA 91125  Department of Physics, California Institute
of Technology, Pasadena, CA 91125  T. Kurumaji  Department of Physics,
Massachusetts Institute of Technology, 77 Massachusetts Ave., Cambridge, MA
02139  G. Refael  Institute for Quantum Information and Matter, California
Institute of Technology, Pasadena, CA 91125  Department of Physics, California
Institute of Technology, Pasadena, CA 91125  J. G. Checkelsky  Department of
Physics, Massachusetts Institute of Technology, 77 Massachusetts Ave.,
Cambridge, MA 02139  D. Hsieh<EMAIL_ADDRESS>Institute for Quantum
Information and Matter, California Institute of Technology, Pasadena, CA 91125
Department of Physics, California Institute of Technology, Pasadena, CA 91125

###### Abstract

In the presence of electron-phonon coupling, an excitonic insulator harbors
two degenerate ground states described by an Ising-type order parameter.
Starting from a microscopic Hamiltonian, we derive the equations of motion for
the Ising order parameter in the phonon coupled excitonic insulator Ta2NiSe5
and show that it can be controllably reversed on ultrashort timescales using
appropriate laser pulse sequences. Using a combination of theory and time-
resolved optical reflectivity measurements, we report evidence of such order
parameter reversal in Ta2NiSe5 based on the anomalous behavior of its
coherently excited order-parameter-coupled phonons. Our work expands the field
of ultrafast order parameter control beyond spin and charge ordered materials.

Exploring new pathways to optically switch Ising-type electronic order
parameters is a major theme of current ultrafast science. In recent years, a
variety of out-of-equilibrium protocols have been developed for rapidly
switching ferromagnetic Lambert  _et al._ (2014); Tudosa  _et al._ (2004);
Gerrits  _et al._ (2002), ferrimagnetic Stanciu  _et al._ (2007); Vahaplar
_et al._ (2009); Radu  _et al._ (2011), antiferromagnetic Kimel  _et al._
(2009); Manz  _et al._ (2016); Schlauderer  _et al._ (2019), and ferroelectric
Fahy and Merlin (1994); Qi  _et al._ (2009); Mankowsky  _et al._ (2017) order
parameters. However, far less is understood about the mechanisms for switching
more exotic order parameters that are not of magnetic and charge dipolar type.

A particularly interesting case is the excitonic insulator (EI), a strongly
correlated electronic phase realized through condensation of bound electron-
hole pairs Jérome  _et al._ (1967). The free energy landscapes of the complex
electronic order parameter and the real lattice order parameter of an EI are
typically characterized by a Mexican hat with continuous U(1) symmetry and a
parabola, respectively [Fig. 1(a)]. However, strong electron-phonon coupling
(EPC) induces a tilting of both the lattice and electronic potentials Lian
_et al._ (2020); Larkin  _et al._ (2017); Kaneko  _et al._ (2013); Nakano  _et
al._ (2018); Suzuki  _et al._ , reducing the U(1) symmetry to a discrete
$Z_{2}$ Ising-type symmetry. Like in magnetic or charge dipole ordered ferroic
materials, this leads to two degenerate ground states characterized by order
parameters of equal magnitude but opposite phase [Fig. 1(b)].

There is currently no experimental method to switch nor directly measure the
phase of an EI order parameter on ultrashort timescales. An alternative
strategy is to measure the phase of the coupled structural order parameter.
However, existing time-resolved x-ray and electron diffraction techniques are
not phase sensitive. Optical phase-resolved second harmonic generation
measurements Mankowsky  _et al._ (2017); Fiebig  _et al._ (2004, 1995) have
been used to measure the phase of structural order parameters in
noncentrosymmetric ferroic materials Fiebig  _et al._ (1994, 2002); Van Aken
_et al._ (2007); Matsubara  _et al._ (2015), but all known EI candidates are
centrosymmetric Kogar  _et al._ (2017); Wakisaka  _et al._ (2009). The
possible presence of 180∘ EI domains further complicates such measurements
because domain averaging would cause overall signal cancellation.

In this Letter, we demonstrate via theory and experiment a pathway to
optically switch an EI order parameter and to probe this reversal through a
coherent EI order-parameter-coupled phonon (OPCP). The time evolution of the
coupled EI and structural order parameters following impulsive laser
excitation are derived from modeling a prototypical system Ta2NiSe5, which we
assume to harbor an EI phase, with an elementary spinless two-band
Hamiltonian. Our simulations reveal that the EI order parameter is stably
reversed above a critical laser fluence, identifiable indirectly via a
saturation of the coherent OPCP amplitude. Using time-resolved coherent phonon
spectroscopy measurements, we experimentally verify this scenario and also
demonstrate how switching can be controlled through the relative timing
between successive laser excitation pulses.

The quasi-one-dimensional (1D) direct band-gap semiconductor Ta2NiSe5 is
reported to undergo an EI transition at a critical temperature $T_{c}=$ 328 K
Wakisaka  _et al._ (2009); Seki  _et al._ (2014), accompanied by a weak
orthorhombic-to-monoclinic structural distortion due to EPC Kaneko  _et al._
(2013). Impulsive laser excitation below $T_{c}$ has been shown to coherently
excite at least five distinct Raman-active phonons with frequencies near 1, 2,
3, 3.7, and 4 THz. The 1, 2, and 4 THz modes are sensitive to the EI
transition at $T_{c}$ Werdehausen  _et al._ (2018); Mor  _et al._ (2018);
Suzuki  _et al._ ; Tang  _et al._ and thus constitute the OPCPs, while the 3
and 3.7 THz modes are reportedly not coupled to the EI order parameter and
thus serve as a control.

Figure 1: Schematic of the electronic and structural free energy landscapes
(a) without and (b) with EPC. In the latter case, pulsed excitation can drive
the system between two degenerate ground states. (c) Schematic of the 1D
spinless two-band model. Red (blue) circles on each site $i$ denote Ta 5$d$
conduction band (Ni 3$d$ valence band) states. The microscopic parameters
discussed in the main text are defined pictorially.

The low energy electronic structure of Ta2NiSe5 consists of a conduction band
with Ta 5$d$ orbital character and a valence band with Ni 3$d$-Se 4$p$
hybridized orbital character. The EI instability is well captured by a 1D
spinless two-band Hamiltonian with EPC Zenker  _et al._ (2014); Tanaka  _et
al._ (2018); Murakami  _et al._ (2017); Tanabe  _et al._ (2018) [Fig. 1(c)],

$\begin{split}H=\sum_{k}(\epsilon_{k}{c_{k}}^{\dagger}c_{k}+\mu_{k}{v_{k}}^{\dagger}v_{k})+\sum_{i}[V{c_{i}}^{\dagger}c_{i}{v_{i}}^{\dagger}v_{i}\\\
+{\omega}_{0}{b_{i}}^{\dagger}b_{i}+g({b_{i}}^{\dagger}+b_{i})({c_{i}}^{\dagger}v_{i}+{v_{i}}^{\dagger}c_{i})],\end{split}$
(1)

where ${c_{k}}^{\dagger},c_{k}$ and ${v_{k}}^{\dagger},v_{k}$ are the
fermionic creation and annihilation operators for conduction and valence band
electrons with momentum $k$, respectively, and ${b_{i}}^{\dagger},b_{i}$ are
the bosonic creation and annihilation operators for an OPCP mode of energy
${\omega}_{0}$ at site $i$. The conduction and valence band dispersions are
given by ${\epsilon}_{k}=\frac{\Delta}{2}+2J_{c}{\sin{(ka/2)}}^{2}$ and
${\mu}_{k}=-\frac{\Delta}{2}-2J_{v}{\sin{(ka/2)}}^{2}$, respectively, where
$2J_{c}$ and $2J_{v}$ are their bandwidths and $\Delta$ is the band gap
Murakami  _et al._ (2017); Tanabe  _et al._ (2018). An on-site interband
electron-electron interaction term $V$ drives the excitonic pairing, and an
EPC term $g$ couples the electronic and phononic subsystems. These microscopic
parameters have been experimentally determined Murakami  _et al._ (2017);
Tanabe  _et al._ (2018).

From Eq.(1), we derive the equations of motion for the EI and structural order
parameters, defined as ${\Phi}_{i}={\langle}{c_{i}}^{\dagger}v_{i}{\rangle}$
and $X_{i}={\langle}{b_{i}}^{\dagger}+b_{i}{\rangle}$, respectively, in two
steps. We first derive the exact expression for the nonequilibrium free energy
functional of $\Phi$ and $X$ in the Keldysh path integral framework, where we
include the light excitation via Peierls substitution. Then we obtain the
equations of motion as the saddle point of the free energy, ignoring
population in the conduction bands, spatial fluctuations of the order
parameter field, and higher-order contributions $O(\Phi^{6})$ (see
Supplemental Material SM ),

$i{\partial}_{t}{\Phi}=[-D(\nabla-
iqA)^{2}+m+U{|{\Phi}|}^{2}]{\Phi}+2g^{\prime}X,$ (2)
${{\partial}_{t}}^{2}{X}=-({{\omega}_{0}}^{2}+\frac{2g^{2}{\omega}_{0}}{V})X-2g{\omega}_{0}\operatorname{Re}(\Phi).$
(3)

Here $D$ is an effective diffusion coefficient, $q$ is the electron charge,
$A$ is the light vector potential, $g^{\prime}$ is the renormalized EPC
coefficient for the electronic channel, and $m$ and $U$ are the second- and
fourth-order expansion coefficients of the electron-electron interaction term,
respectively. These parameters are all functions of $J_{c}$, $J_{v}$,
$\Delta$, $g$, and $V$ SM . We further introduce two phenomenological
constants $\gamma_{e}$ and $\gamma_{ph}$ to account for damping of electronic
and structural modes due to thermal fluctuation, whose values can be
experimentally determined SM .

Figure 2: Simulation results of $V(\Phi)$ and $V(X)$ (defined in main text)
for $F>F_{c}$ with experimentally determined parameters SM . Snapshots of the
potential landscapes (solid lines) and the electronic and structural order
parameters (circles) are shown (a) in the equilibrium state (reproduced as
dashed lines in (b)-(d)), (b) at the moment of excitation, (c) during transit
into the reversed state, and (d) in the reversed state before equilibration,
where both potentials are modulated at the phonon frequency (red arrows).
Axes’ scales are the same in all panels.

We first present a qualitative picture of how order parameter reversal occurs
in our model. In equilibrium the electronic potential
$V(\Phi)=\frac{1}{2}m{\Phi}^{2}+\frac{1}{4}U{\Phi}^{4}+2g^{\prime}X{\Phi}$ has
either a tilted parabolic ($m>$ 0) or tilted Mexican-hat ($m<$ 0) form
capturing the absence or presence of exciton condensation. For homogeneous
optical excitation, one can ignore spatial derivatives of $\Phi$. For pulsed
excitation with infrared light, whose frequency well exceeds the electronic
Higgs-Goldstone Murakami  _et al._ (2017, 2020) and OPCP mode frequencies, one
can also average out the fast oscillations of the perturbation and retain only
its Gaussian envelope. Under these conditions one can make the simplification
$m(t)=-D(\nabla-iqA)^{2}+m$ ${\triangleq}f(t)+m$, where
$f(t)={\alpha}F\exp{(-\frac{4\ln({2})t^{2}}{\sigma^{2}})}$. Here $\sigma$ is
the temporal width of the Gaussian pulse, $\alpha$ is a positive constant
scaling factor that can be calculated analytically SM , and $F$ is the pump
fluence, the only tunable parameter in our model. In the EI phase ($m<$ 0),
optical excitation therefore acts to instantaneously increase $m(t)$. The
subsequent reduction of the EI order parameter, which occurs on a timescale
much shorter than $2\pi/\omega_{0}$, results in a sudden shift in the lattice
potential
$V(X)=\frac{1}{2}{\omega_{0}}^{2}{X}^{2}+2g{\omega_{0}}\operatorname{Re}(\Phi)X$
due to EPC SM , launching coherent oscillations through displacive excitation.
In the low fluence regime, where the phonon oscillation amplitude is small
enough such that $X$ does not change sign, the direction of tilt of both
$V(\Phi)$ and $V(X)$ remains unchanged and so no switching occurs. However,
above a critical fluence $F_{c}$, where the phonon oscillation amplitude
becomes large enough to change the sign of $X$, the tilting of both potentials
is reversed and the system can relax into the switched state.

Numerical simulations of our model using experimentally determined material
parameters for Ta2NiSe5 and $\sigma$ = 100 fs were carried out in the EI phase
SM . Figure 2 displays simulation results for $F$ slightly greater than
$F_{c}$. At the moment of excitation $t=0$ [Fig. 2(b)], there is an
instantaneous change in $V(\Phi)$ from Mexican-hat to parabolic form caused by
the light-induced enhancement of $m(t)$. The EI order parameter evolves
rapidly to the new potential minimum with overdamped dynamics and is quenched
within the pulse duration. This leads to a rapid shift in $V(X)$, launching
coherent oscillations of the underdamped OPCP, which shakes the electronic
potential via EPC at the phonon frequency. Once the pulsed excitation is over,
$V(\Phi)$ recovers a Mexican-hat form. However, as $X$ crosses zero within the
first half period of oscillation [Fig. 2(c)], the tilting of $V(\Phi)$ is
reversed with respect to the pre-pumped ($t<$ 0) case [Fig. 2(a)], sending
$\Phi$ toward the new minimum on the negative side. As $\Phi$ crosses zero,
the tilt of $V(X)$ is also reversed due to EPC, thus pushing $X$ to the new
minimum on the positive side. The system then continues to oscillate about the
reversed minima at the phonon frequency until the OPCP is damped out [Fig.
2(d)]. Note that this model treats light as a coherent drive without
considering heating- and cooling-induced changes in $m(t)$. However,
accounting for the latter merely shifts $F_{c}$ SM .

In lieu of probing the phase of $X$ and $\Phi$, we propose that the reversal
can be identified via the pump fluence dependence of the OPCP amplitude. We
solve Eqs. (2) and (3) to obtain the time evolution of $\Phi$ and $X$ and then
define the OPCP amplitude by its peak height in the fast Fourier transform
(FFT) of $X$. A conventional Raman-active phonon is coherently launched
through either displacive excitation or impulsive stimulated Raman scattering
(ISRS) Stevens  _et al._ (2002); Melnikov  _et al._ (2011) with an amplitude
that is linearly proportional to $F$. For $F<F_{c}$, the amplitude of an OPCP
also scales linearly with $F$. In this regime the structural order parameter
simply oscillates about the initial potential minimum and thus behaves like a
conventional phonon. But once the initial displacement of $V(X)$ is large
enough to enable escape to the opposite minimum ($F>F_{c}$), the amplitude
ceases to grow. Order parameter reversal is thus marked by a saturation in the
amplitude versus fluence curve.

We argue that experimental evidence for this phenomenon already exists in
published studies of Ta2NiSe5. Werdehausen $et$ $al.$ Werdehausen  _et al._
(2018) performed ultrafast optical reflectivity measurements at $T$ = 120 K
using a pump photon energy of 1.55 eV and observed clear coherent oscillations
of the 1 THz OPCP as well as the uncoupled 3 THz phonon (ISRS). The pump
fluence dependence of these two mode amplitudes is reproduced in Fig. 3. The 3
THz mode scales linearly with $F$, consistent with its assignment as a
conventional phonon. In contrast, the 1 THz mode scales linearly with $F$ only
at low fluences and then saturates above $\sim$0.4 mJ/cm2, consistent with an
order parameter reversal. By overlaying our simulation results SM  atop these
experimental curves, we find close agreement (Fig. 3).

Our theory also predicts that the exciton condensate should be transiently
quenched [$m(t)\rightarrow 0$] above a critical fluence $F=F^{*}$ where the
condition $f(t)=|m|$ is satisfied. Eqs. (2) and (3) do not constrain $F^{*}$
to coincide with $F_{c}$ and our simulation shows that $F^{*}$ is clearly
lower than $F_{c}$ in Ta2NiSe5 (Fig. 3). Recently Tang $et$ $al.$ Tang  _et
al._ performed time- and angle-resolved photoemission spectroscopy (tr-ARPES)
measurements on Ta2NiSe5 at $T$ = 30 K and tracked the dynamics of the charge
gap, a measure of $\Phi$, immediately after pumping with 1.77 eV light
polarized perpendicular to the chain direction (equivalent geometry to Ref.
Werdehausen  _et al._ (2018)). They found that the instantaneous gap size
decreases linearly with increasing pump fluence and saturates above 0.29
mJ/cm2, which was interpreted as the point where $\Phi$ transiently collapses.
The fact that this fluence is lower than 0.4 mJ/cm2, and is expected to be
even lower if the experiment were conducted at 120 K, is consistent with our
theory.

Figure 3: Experimental pump fluence dependence of the 1 THz (red circles) and
3 THz (blue circles) coherent phonon amplitudes in Ta2NiSe5 reproduced from
Ref.Werdehausen  _et al._ (2018). Simulation results for an OPCP (red line)
and a conventional ISRS phonon (blue line) are overlaid and horizontally
scaled ($\alpha\approx 700$) to match the experimental data. Vertical dashed
lines mark the calculated $F_{c}$ and $F^{*}$. The non-monotonic behavior of
the OPCP amplitude just above $F_{c}$ arises from strong feedback between
$\Phi$ and $X$ immediately after excitation SM .

The dynamical nature of the order parameter reversal process suggests that it
can be controlled not merely by the total pump energy deposited, but also by
its distribution in time. To show this, we consider a situation where the
sample is pumped by two identical pulses separated by time $\delta t$, with
individual fluences $F<F_{c}$ but $2F>F_{c}$. For $\delta t\rightarrow 0$ the
system is effectively pumped by a single pulse exceeding $F_{c}$ and so
reversal occurs, while for $\delta t\rightarrow\infty$, the system relaxes
back to the initial ground state before the second pulse arrives and so no
reversal occurs. To qualitatively understand the behavior at intermediate
$\delta t$ values where the system is still dynamically evolving when the
second pulse arrives, we recall that in the single pulse case, switching
occurs once the OPCP amplitude is large enough to change the sign of $X$.
Therefore, in the two-pulse case, switching possibly occurs if coherent
oscillations of the OPCP induced by the first pulse can be sufficiently
amplified by the second time-delayed pulse.

Figure 4: Simulated time evolution of (a) Re($\Phi$) and (b) $X$ following
single-pulse pumping (black) and two-pulse OP (blue) and IP (red) pumping of
the 2 THz phonon using $F=0.96F_{c}$ and the same microscopic parameters as in
Figs. 2 and 3. The calculated instantaneous electronic potential is displayed
at several select times. (c) Normalized FFT of the traces shown in (b). Each
curve is normalized by the peak value of the single-pulse pumping curve. (d)
Reflectivity transients measured from Ta2NiSe5 under the same pumping
conditions used in the simulations. Curves are vertically offset for clarity.
(e) Normalized FFT of the traces shown in (d).

Previous studies have shown that impulsively and displacively excited
conventional Raman-active phonons can be coherently amplified (suppressed) by
a second pump pulse when $\delta t$ is an integer (half-integer) multiple of
the phonon period Hase  _et al._ (1996, 1998); Murray  _et al._ (2005);
Misochko  _et al._ (2007); Wu and Xu (2007); Cheng  _et al._ (2017), dubbed
in-phase (IP) and out-of-phase (OP) pumping respectively. Therefore we
simulated the effects of both IP and OP pumping on the order parameters of
Ta2NiSe5 using the same material parameters as before. We chose to simulate
the OPCP at $\omega_{0}/2\pi$ = 2 THz rather than at 1 THz because recent tr-
ARPES data show that the most pronounced modulations of the valence band
maximum occur at 2 THz, suggesting strong coupling to $\Phi$ Suzuki  _et al._
. As shown in Fig. 4(a), pumping by a single pulse with $F$ slightly less than
$F_{c}$ causes a rapid but incomplete reduction of $\Phi$, followed by a
slower recovery back to its original value on a timescale set by the damping
of the 2 THz phonon. For the two-pulse case, OP pumping of the 2 THz phonon
similarly leads to partial suppression of $\Phi$ without reversal, but
reversal is achieved with IP pumping. This phenomenon is again manifested
through an unconventional behavior of the OPCP. As shown in Figs. 4(b) and
(c), OP pumping leads to suppression of the 2 THz phonon amplitude relative to
the single pump case, resembling a conventional phonon because the oscillation
is around the initial potential minimum. But, in contrast to conventional
behavior, IP pumping does not lead to further amplification once $X$ is
excited to the opposite minimum.

To verify this prediction, we performed transient optical reflectivity
measurements on Ta2NiSe5 single crystals SM  using two identical pump pulses
($\sigma$ = 80 fs) with variable $\delta t$. The light was polarized
perpendicular to the chain direction and the fluence of each pulse was tuned
slightly below $F_{c}$ to match our simulations. We chose a pump photon energy
of 1 eV to enhance the 2 THz oscillations SM . Figure 4(d) shows the
fractional change in reflectivity ($\Delta R/R$) versus time for both IP and
OP pumping of the 2 THz phonon, as well as for pumping with only a single
pulse. All three curves exhibit fast ($\sim$1 ps) exponential decay following
pump excitation, corresponding to the charge relaxation process. Oscillations
from the beating of several coherently excited phonons are also clear. A FFT
of the data shows the most pronounced peaks at 2, 3, and 3.8 THz SM . A focus
on the 2 THz mode reveals that OP pumping strongly suppresses its amplitude
relative to the single-pump case whereas IP pumping does not amplify it [Fig.
4(e)], in quantitative agreement with our simulations [Fig. 4(c)]. In
contrast, strong amplification occurs for the 3 and 3.8 THz modes SM ,
consistent with their uncoupled nature.

Our field theory description of the EI order parameter goes beyond the
phenomenological time-dependent Landau theory Huber  _et al._ (2014); Zong
_et al._ (2019); Yusupov  _et al._ (2010); Harter  _et al._ (2018) in that it
allows the order parameter to explore the tilted Mexican-hat potential in the
complex plane and can be naturally linked to microscopic parameters of the
underlying lattice model. While more details including extension beyond the
mean-field limit ($\nabla\neq 0$), temperature dependence with $T>0$,
diffusion perpendicular to the surface Yusupov  _et al._ (2010), time-
dependent damping Huber  _et al._ (2014), and anharmonic phonon coupling Hase
_et al._ (2002) can be added to refine the simulations, our minimal
microscopic theory already captures the most salient physics and experimental
features. These ideas and dynamical protocols apply not only to excitonic
insulators, but also to any system featuring a continuous-symmetry-breaking
electronic order parameter induced by coupling to a structural order
parameter, such as a charge ordered system coupled to a Peierls distortion or
an orbital ordered system coupled to a Jahn-Teller distortion. Therefore the
OPCP behavior revealed here may be a general fingerprint of electronic order
parameter switching.

###### Acknowledgements.

We thank Rick Averitt, Edoardo Baldini, Swati Chaudhary, Nuh Gedik, Xinwei Li,
Tianwei Tang and Alfred Zong for useful discussions. D.H. and J.G.C.
acknowledge support from the DARPA DRINQS program (Grant No. D18AC00014). D.H.
also acknowledges support for instrumentation from the Institute for Quantum
Information and Matter, a NSF Physics Frontiers Center (PHY-1733907). M.B.
acknowledges the support from the Department of Energy under Award No. DE-
SC0019166. T.K. acknowledges the support by the Yamada Science Foundation
Fellowship for Research Abroad and JSPS Overseas Research Fellowships.

## References

  * Lambert  _et al._ (2014) C.-H. Lambert, S. Mangin, B. S. D. C. S. Varaprasad, Y. K. Takahashi, M. Hehn, M. Cinchetti, G. Malinowski, K. Hono, Y. Fainman, M. Aeschlimann, and E. E. Fullerton, Science 345, 1337 (2014).
  * Tudosa  _et al._ (2004) I. Tudosa, C. Stamm, A. B. Kashuba, F. King, H. C. Siegmann, J. Stöhr, G. Ju, B. Lu, and D. Weller, Nature 428, 831 (2004).
  * Gerrits  _et al._ (2002) T. Gerrits, H. A. M. van den Berg, J. Hohlfeld, L. Bär, and T. Rasing, Nature 418, 509 (2002).
  * Stanciu  _et al._ (2007) C. D. Stanciu, F. Hansteen, A. V. Kimel, A. Kirilyuk, A. Tsukamoto, A. Itoh, and T. Rasing, Phys. Rev. Lett. 99, 047601 (2007).
  * Vahaplar  _et al._ (2009) K. Vahaplar, A. M. Kalashnikova, A. V. Kimel, D. Hinzke, U. Nowak, R. Chantrell, A. Tsukamoto, A. Itoh, A. Kirilyuk, and T. Rasing, Phys. Rev. Lett. 103, 117201 (2009).
  * Radu  _et al._ (2011) I. Radu, K. Vahaplar, C. Stamm, T. Kachel, N. Pontius, H. A. Dürr, T. A. Ostler, J. Barker, R. F. L. Evans, R. W. Chantrell, A. Tsukamoto, A. Itoh, A. Kirilyuk, T. Rasing, and A. V. Kimel, Nature 472, 205 (2011).
  * Kimel  _et al._ (2009) A. V. Kimel, B. A. Ivanov, R. V. Pisarev, P. A. Usachev, A. Kirilyuk, and T. Rasing, Nature Physics 5, 727 (2009).
  * Manz  _et al._ (2016) S. Manz, M. Matsubara, T. Lottermoser, J. Büchi, A. Iyama, T. Kimura, D. Meier, and M. Fiebig, Nature Photonics 10, 653 (2016).
  * Schlauderer  _et al._ (2019) S. Schlauderer, C. Lange, S. Baierl, T. Ebnet, C. P. Schmid, D. C. Valovcin, A. K. Zvezdin, A. V. Kimel, R. V. Mikhaylovskiy, and R. Huber, Nature 569, 383 (2019).
  * Fahy and Merlin (1994) S. Fahy and R. Merlin, Phys. Rev. Lett. 73, 1122 (1994).
  * Qi  _et al._ (2009) T. Qi, Y.-H. Shin, K.-L. Yeh, K. A. Nelson, and A. M. Rappe, Phys. Rev. Lett. 102, 247603 (2009).
  * Mankowsky  _et al._ (2017) R. Mankowsky, A. von Hoegen, M. Först, and A. Cavalleri, Phys. Rev. Lett. 118, 197601 (2017).
  * Jérome  _et al._ (1967) D. Jérome, T. M. Rice, and W. Kohn, Phys. Rev. 158, 462 (1967).
  * Lian  _et al._ (2020) C. Lian, S.-J. Zhang, S.-Q. Hu, M.-X. Guan, and S. Meng, Nature Communications 11, 43 (2020).
  * Larkin  _et al._ (2017) T. I. Larkin, A. N. Yaresko, D. Pröpper, K. A. Kikoin, Y. F. Lu, T. Takayama, Y.-L. Mathis, A. W. Rost, H. Takagi, B. Keimer, and A. V. Boris, Phys. Rev. B 95, 195144 (2017).
  * Kaneko  _et al._ (2013) T. Kaneko, T. Toriyama, T. Konishi, and Y. Ohta, Phys. Rev. B 87, 035121 (2013).
  * Nakano  _et al._ (2018) A. Nakano, T. Hasegawa, S. Tamura, N. Katayama, S. Tsutsui, and H. Sawa, Phys. Rev. B 98, 045139 (2018).
  * (18) T. Suzuki, Y. Shinohara, Y. Lu, M. Watanabe, J. Xu, K. L. Ishikawa, H. Takagi, M. Nohara, N. Katayama, H. Sawa, M. Fujisawa, T. Kanai, J. Itatani, T. Mizokawa, S. Shin, and K. Okazaki, arXiv:2002.10037 .
  * Fiebig  _et al._ (2004) M. Fiebig, T. Lottermoser, D. Fröhlich, and S. Kallenbach, Opt. Lett. 29, 41 (2004).
  * Fiebig  _et al._ (1995) M. Fiebig, D. Fröhlich, G. Sluyterman v. L., and R. V. Pisarev, Applied Physics Letters 66, 2906 (1995).
  * Fiebig  _et al._ (1994) M. Fiebig, D. Fröhlich, B. B. Krichevtsov, and R. V. Pisarev, Phys. Rev. Lett. 73, 2127 (1994).
  * Fiebig  _et al._ (2002) M. Fiebig, T. Lottermoser, D. Fröhlich, A. V. Goltsev, and R. V. Pisarev, Nature 419, 818 (2002).
  * Van Aken  _et al._ (2007) B. B. Van Aken, J.-P. Rivera, H. Schmid, and M. Fiebig, Nature 449, 702 (2007).
  * Matsubara  _et al._ (2015) M. Matsubara, S. Manz, M. Mochizuki, T. Kubacka, A. Iyama, N. Aliouane, T. Kimura, S. L. Johnson, D. Meier, and M. Fiebig, Science 348, 1112 (2015).
  * Kogar  _et al._ (2017) A. Kogar, M. S. Rak, S. Vig, A. A. Husain, F. Flicker, Y. I. Joe, L. Venema, G. J. MacDougall, T. C. Chiang, E. Fradkin, J. van Wezel, and P. Abbamonte, Science 358, 1314 (2017).
  * Wakisaka  _et al._ (2009) Y. Wakisaka, T. Sudayama, K. Takubo, T. Mizokawa, M. Arita, H. Namatame, M. Taniguchi, N. Katayama, M. Nohara, and H. Takagi, Phys. Rev. Lett. 103, 026402 (2009).
  * Seki  _et al._ (2014) K. Seki, Y. Wakisaka, T. Kaneko, T. Toriyama, T. Konishi, T. Sudayama, N. L. Saini, M. Arita, H. Namatame, M. Taniguchi, N. Katayama, M. Nohara, H. Takagi, T. Mizokawa, and Y. Ohta, Phys. Rev. B 90, 155116 (2014).
  * Werdehausen  _et al._ (2018) D. Werdehausen, T. Takayama, M. Höppner, G. Albrecht, A. W. Rost, Y. Lu, D. Manske, H. Takagi, and S. Kaiser, Sci. Adv. 4 (2018), eaap8652.
  * Mor  _et al._ (2018) S. Mor, M. Herzog, J. Noack, N. Katayama, M. Nohara, H. Takagi, A. Trunschke, T. Mizokawa, C. Monney, and J. Stähler, Phys. Rev. B 97, 115154 (2018).
  * (30) T. Tang, H. Wang, S. Duan, Y. Yang, C. Huang, Y. Guo, D. Qian, and W. Zhang, arXiv:2003.00514 .
  * Zenker  _et al._ (2014) B. Zenker, H. Fehske, and H. Beck, Phys. Rev. B 90, 195118 (2014).
  * Tanaka  _et al._ (2018) Y. Tanaka, M. Daira, and K. Yonemitsu, Phys. Rev. B 97, 115105 (2018).
  * Murakami  _et al._ (2017) Y. Murakami, D. Golež, M. Eckstein, and P. Werner, Phys. Rev. Lett. 119, 247601 (2017).
  * Tanabe  _et al._ (2018) T. Tanabe, K. Sugimoto, and Y. Ohta, Phys. Rev. B 98, 235127 (2018).
  * (35) See Supplemental Material for extensive derivation of theory and simulation and experimental details, which contains Refs. [36-45].
  * Kamenev (2011) A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, 2011).
  * Mor  _et al._ (2017) S. Mor, M. Herzog, D. Golež, P. Werner, M. Eckstein, N. Katayama, M. Nohara, H. Takagi, T. Mizokawa, C. Monney, and J. Stähler, Phys. Rev. Lett. 119, 086401 (2017).
  * Okazaki  _et al._ (2018) K. Okazaki, Y. Ogawa, T. Suzuki, T. Yamamoto, T. Someya, S. Michimae, M. Watanabe, Y. Lu, M. Nohara, H. Takagi, N. Katayama, H. Sawa, M. Fujisawa, T. Kanai, N. Ishii, J. Itatani, T. Mizokawa, and S. Shin, Nature Communications 9, 4322 (2018).
  * (39) M. D. Watson, I. Marković, E. A. Morales, P. L. Fèvre, M. Merz, A. A. Haghighirad, and P. D. C. King, arXiv:1912.01591 .
  * Lu  _et al._ (2017) Y. F. Lu, H. Kono, T. I. Larkin, A. W. Rost, T. Takayama, A. V. Boris, B. Keimer, and H. Takagi, Nature Communications 8, 14408 (2017).
  * Zeiger  _et al._ (1992) H. J. Zeiger, J. Vidal, T. K. Cheng, E. P. Ippen, G. Dresselhaus, and M. S. Dresselhaus, Phys. Rev. B 45, 768 (1992).
  * Kuznetsov and Stanton (1994) A. V. Kuznetsov and C. J. Stanton, Phys. Rev. Lett. 73, 3243 (1994).
  * Misochko and Lebedev (2016) O. V. Misochko and M. V. Lebedev, Phys. Rev. B 94, 184307 (2016).
  * Hunsche  _et al._ (1995) S. Hunsche, K. Wienecke, T. Dekorsy, and H. Kurz, Phys. Rev. Lett. 75, 1815 (1995).
  * Tangney and Fahy (1999) P. Tangney and S. Fahy, Phys. Rev. Lett. 82, 4340 (1999).
  * Murakami  _et al._ (2020) Y. Murakami, D. Golež, T. Kaneko, A. Koga, A. J. Millis, and P. Werner, Phys. Rev. B 101, 195118 (2020).
  * Stevens  _et al._ (2002) T. E. Stevens, J. Kuhl, and R. Merlin, Phys. Rev. B 65, 144304 (2002).
  * Melnikov  _et al._ (2011) A. Melnikov, O. Misochko, and S. Chekalin, Physics Letters A 375, 2017 (2011).
  * Hase  _et al._ (1996) M. Hase, K. Mizoguchi, H. Harima, S. Nakashima, M. Tani, K. Sakai, and M. Hangyo, Applied Physics Letters 69, 2474 (1996).
  * Hase  _et al._ (1998) M. Hase, T. Itano, K. Mizoguchi, and S. ichi Nakashima, Japanese Journal of Applied Physics 37, L281 (1998).
  * Murray  _et al._ (2005) E. D. Murray, D. M. Fritz, J. K. Wahlstrand, S. Fahy, and D. A. Reis, Phys. Rev. B 72, 060301(R) (2005).
  * Misochko  _et al._ (2007) O. V. Misochko, M. V. Lebedev, H. SchÃ€fer, and T. Dekorsy, Journal of Physics: Condensed Matter 19, 406220 (2007).
  * Wu and Xu (2007) A. Q. Wu and X. Xu, Applied Physics Letters 90, 251111 (2007).
  * Cheng  _et al._ (2017) Y.-H. Cheng, F. Y. Gao, S. W. Teitelbaum, and K. A. Nelson, Phys. Rev. B 96, 134302 (2017).
  * Huber  _et al._ (2014) T. Huber, S. O. Mariager, A. Ferrer, H. Schäfer, J. A. Johnson, S. Grübel, A. Lübcke, L. Huber, T. Kubacka, C. Dornes, C. Laulhe, S. Ravy, G. Ingold, P. Beaud, J. Demsar, and S. L. Johnson, Phys. Rev. Lett. 113, 026401 (2014).
  * Zong  _et al._ (2019) A. Zong, P. E. Dolgirev, A. Kogar, E. Ergecen, M. B. Yilmaz, Y.-Q. Bie, T. Rohwer, I.-Cheng Tung, J. Straquadine, X. Wang, Y. Yang, X. Shen, R. Li, J. Yang, S. Park, M. C. Hoffmann, B. K. Ofori-Okai, M. E. Kozina, H. Wen, X. Wang, I. R. Fisher, P. Jarillo-Herrero, and N. Gedik, Phys. Rev. Lett. 123, 097601 (2019).
  * Yusupov  _et al._ (2010) R. Yusupov, T. Mertelj, V. V. Kabanov, S. Brazovskii, P. Kusar, J.-H. Chu, I. R. Fisher, and D. Mihailovic, Nature Physics 6, 681 (2010).
  * Harter  _et al._ (2018) J. W. Harter, D. M. Kennes, H. Chu, A. de la Torre, Z. Y. Zhao, J.-Q. Yan, D. G. Mandrus, A. J. Millis, and D. Hsieh, Phys. Rev. Lett. 120, 047601 (2018).
  * Hase  _et al._ (2002) M. Hase, M. Kitajima, S.-i. Nakashima, and K. Mizoguchi, Phys. Rev. Lett. 88, 067401 (2002).


lemmatheorem lemma corollarytheorem corollary definitiontheorem definition

# Maximizing the Sum of Radii of Disjoint Balls or Disks

David Eppstein Computer Science Department, University of California, Irvine;
Irvine, CA, 92697, USA<EMAIL_ADDRESS>Research supported in part by NSF
grants CCF-1228639, CCF-1618301, and CCF-1616248.

###### Abstract

Finding nonoverlapping balls with given centers in any metric space,
maximizing the sum of radii of the balls, can be expressed as a linear
program. Its dual linear program expresses the problem of finding a minimum-
weight set of cycles (allowing 2-cycles) covering all vertices in a complete
geometric graph. For points in a Euclidean space of any finite dimension $d$,
with any convex distance function on this space, this graph can be replaced by
a sparse subgraph obeying a separator theorem. This graph structure leads to
an algorithm for finding the optimum set of balls in time $O(n^{2-1/d})$,
improving the $O(n^{3})$ time of a naive cycle cover algorithm. As a
subroutine, we provide an algorithm for weighted bipartite matching in graphs
with separators, which speeds up the best previous algorithm for this problem
on planar bipartite graphs from $O(n^{3/2}\log n)$ to $O(n^{3/2})$ time. We
also show how to constrain the balls to all have radius at least a given
threshold value, and how to apply our radius-sum optimization algorithms to
the problem of embedding a finite metric space into a star metric minimizing
the average distance to the hub.

##  1 Introduction

In this paper we consider the following problem. We are given a set of points
in the plane, and must choose a radius for each point, so that the disks with
the given radii do not overlap and the sum of radii is as large as possible.
The same problem can be generalized to arbitrary metric spaces, with some care
about definitions to allow us to work with finite metric spaces (such as a
space whose only points are the given centers). Although we believe these
problems to be natural and interesting on their own merits, we provide two
motivating applications, in map labeling and in optimal embedding of metric
spaces into star metrics.

To formally define the problems we study, we use the following definition.

######  Definition .

We define a _ball_ , in a given metric space, to be a pair $(c,r)$ where $c$
is a point in the metric space (the center of the ball), and $r$ is a non-
negative real number (the radius of the ball). A ball $(c,r)$ contains a point
$p$ if the distance from $c$ to $p$ is less than $r$. Two balls in an
arbitrary metric space _overlap_ when their sum of radii exceeds the center
distance, and are _nonoverlapping_ otherwise. If the sum of radii equals the
center distance, we say that the two balls _touch_.

For a _geodesic_ metric space such as a Euclidean space, in which every two
points can be joined by an isometrically embedded line segment, this
definition coincides with the usual intuition: balls are distinct when they
contain distinct sets of points, and overlap when they have a point of
intersection. However, these definitions allow us to work more conveniently
with non-geodesic metric spaces, such as a finite metric space consisting only
of the ball centers. In such spaces, distinct balls may contain the same sets
of points, and overlapping balls may not have any points that they both
contain.

Then in any metric space, finding a set of nonoverlapping metric balls with
$n$ given center points $p_{i}$, maximizing the sum of non-negative radii
$r_{i}$, can be expressed as a linear program. The objective is to maximize
the linear function $\sum r_{i}$, subject to constraints that each pair of
balls remain nonoverlapping, i.e. that $r_{i}+r_{j}\leq d(p_{i},p_{j})$. This
linear program has two variables per constraint, a well-studied special case
of linear programming. But although strongly polynomial algorithms for
feasibility with two variables per constraint are known [4, 14], they do not
extend to optimization, and their running time is higher than might be
desired. Therefore, it remains of interest to find a purely combinatorial
algorithm for the problem, with as low a running time as possible.

In this paper, we provide such a combinatorial algorithm, running in cubic
time for general metrics and subquadratic time for low-dimensional Euclidean
spaces.

Figure 1: The moat-growing method for Euclidean matching (left) versus the set
of disks with maximum sum of radii centered at the same set of points (right).

The problem that we study here should be distinguished from the classical
“moat-growing” method for solving Euclidean minimum weight matching problems,
in which we grow disks around points until they touch but then, around any
subset of an odd number of points with touching regions, continue growing a
“moat”. This method was illustrated, for instance, on the cover of a
combinatorial optimization book by Cook et al. [5]. In contrast, our method
has only disks, with no moats, and when it finds a subsystem of an odd number
of points with touching disks, those disks are frozen in place, with no
possibility of future growth (section 5). Although similar in visual
appearance to moat-growing (Figure 1), the solution to our problem can have a
substantially smaller value than the solution to the Euclidean minimum weight
matching problem. The radius-sum optimization problem that we consider here
was formulated earlier by Jünger and Pulleyblank [18], as a type of geometric
dual problem to Euclidean minimum weight matching, but they quickly dismissed
it as inaccurate before formulating tighter moat-growing-based heuristics for
matching, and they did not provide an algorithm for finding the optimal radii.

###  1.1 New results

We prove the following results:

  * •

Finding metric balls with maximum sum of radii is equivalent under linear
programming duality to finding a minimum-length set of cycles (allowing
$2$-cycles) that cover all vertices of the complete geometric graph on the
given centers. The maximum sum of radii equals half of the minimum total cycle
length. By reducing cycle covers to weighted matching, the optimal set of
balls in any metric space can be constructed in cubic time.

  * •

For points in Euclidean spaces of bounded dimension $d$, the edges of the
optimal cycle cover can be found in a subgraph of the complete geometric
graph, the intersection graph of nearest neighbor balls of the points. This
graph is sparse, having $O(n)$ edges with a constant of proportionality
depending singly-exponentially in the dimension. Moreover, it obeys a
separator theorem with separators of size $O(n^{1-1/d})$, and the graph and
its separator decomposition can be constructed in time $O(n\log n)$.

  * •

A separator-based divide and conquer weighted matching algorithm can find an
optimal cycle cover and set of nonoverlapping balls in time $O(n^{2-1/d})$. In
particular, we can find an optimal set of disks in the plane in time
$O(n^{3/2})$.

  * •

As a subroutine, we use a new weighted bipartite matching algorithm for
families of graphs that obey a separator theorem, which runs in time
$O(n^{3/2})$ for planar bipartite graphs, improving a previous bound of
$O(n^{3/2}\log n)$ for this case.

  * •

We show how to solve radius-sum optimization problems in which all disks or
balls have a minimum radius greater than zero, by transforming these problems
into unconstrained radius-sum optimization problems on a different metric
space.

To avoid issues of numerical precision we consider only strongly-polynomial-
time algorithms, in a model of computation in which distance computation,
arithmetic, and comparisons take constant time per operation.

###  1.2 Related research

Along with the map labeling research discussed below, researchers have studied
other types of geometric optimization problems in which the optimization
criterion is a sum of radii of balls or circles. These include finding a set
of $k$ balls with given centers, drawn from a larger set of $n$ points, that
cover all points and minimize the sum of radii of the balls [12], finding a
connected set of disks in the plane with given centers that minimize the sum
of radii [2], and finding both a collection of disks centered at a subset of
input points that covers all input points, and a tour connecting the disk
centers, minimizing a combination of disk radii and tour length [1].

##  2 Applications

###  2.1 Map labeling

_Map labeling_ is an algorithmic problem in which one must place non-
overlapping text labels on maps or other visualizations. Researchers of map
labeling have studied various problems of assigning shapes of maximum size to
a given set of points in the plane [6, 26, 17, 21]. One simple case of this
problem, finding circles centered at the given points that maximize the
minimum radius, can be solved by finding the closest pair of points and
setting all radii equal to half of this pair’s distance. However, this measure
of solution quality penalizes the label sizes even for points far away from
the closest pair, where larger radii could be used without overlap (Figure 2,
left). On the other hand, an $L_{2}$ measure of solution quality (maximizing
the sum of disk areas) would also be unsatisfactory: even with only two points
to be labeled, the $L_{2}$ solution would assign zero radius to one point
(Figure 2, center). The problem that we study in this paper is the $L_{1}$
version of this problem, in which we are given as input $n$ points and must
maximize the sum of radii of disjoint disks centered at those points (Figure
2, right). Two points can be optimally labeled by two equal-radius disks, and
this criterion has the largest value of $p$ among $L_{p}$ criteria (maximizing
sums of $p$th powers of the radii) that allow this equal-radius solution.

Figure 2: For disks with fixed centers, maximizing the minimum radius (left)
leads to disks that are smaller than necessary, while maximizing the sum of
areas (center) leads to degenerate zero-radius disks. Maximizing the sum of
radii (right) avoids both problems.

We remark, however, that the algorithms we present here do not provide a
complete solution to the map labeling problem. One reason for this is that
there may be multiple optimal families of disks for a given collection of
points, all with the same sum of radii, and that map labeling would require us
to select a single radius for each point. For instance, with only two points,
any pair of radii whose sum os the distance between the points would be
optimal; in this case, the preferred solution for map labeling would likely be
the one in which the two radii are equal. We leave the problem of how to
formulate the best choice among the optimum solutions to our radius sum
maximization problem, and of finding this best choice by a fast combinatorial
algorithm, as open for future research.

Another, more serious, problem with applying radius sum optimization to map
labeling is that, for more than two disks, the optimal solution may still
involve disks of zero radius. In particular, when the input consists of three
collinear points, the optimal solution assigns zero radius to the middle
point, with the other two points having radii equal to their distance from the
middle point (Figure 3). One way to work around this issue, and ensure that
all disks have non-zero radius, is to constrain all of the disks to have
radius at least $\delta$, for some threshold value $\delta>0$. We may choose
$\delta$ to be as large as half of the minimum distance between the given
points. We then seek a system of radii $r_{i}\geq\delta$ that form non-
overlapping disks of maximum total radius. As we show in section 7, this
constrained version of the problem, in which we have an additional lower bound
on the radii of the balls, can be solved efficiently by transforming it into
an unconstrained problem on a different distance metric.

Figure 3: The maximum radius sum for three collinear points assigns zero
radius to the middle point.

###  2.2 Metric embedding into stars

There has been extensive research on embedding complex metrics into simpler
metrics with low distortion; such methods have many applications in
approximation algorithms, by allowing approximations designed for the simpler
metric to be applied to the more complex one [15]. Eppstein and Wortman [10]
consider one such problem: given any finite metric space $(X,d)$, find an
embedding of it into a star metric space with minimum distortion. Here, a star
is a space in which there is a central hub point $h$, not necessarily one of
the given points, such that the distance between every two points is the sum
of their distances to the hub. The space of all distance-minimal non-
contractive mappings from a given metric space $(X,d)$ to star metrics is
called the _tight span_ of $(X,d)$. With the sup-norm (the maximum difference
between hub distances) it is itself a metric space, and includes an isometric
copy of the original space $(X,d)$ (where each $x$ in $X$ is mapped to the
star having $x$ as a hub) [16]. The problem of Eppstein and Wortman is to
select the optimal point of the tight span, measuring the quality of each
point by the distortion of the corresponding star metric. We note that this is
not the same as choosing an optimal hub point from some ambient space
containing the input; for instance, if the input is the four points of a unit
square in the Euclidean plane, the optimal star metric gives each of the four
input points distance $1/2$ from the hub, while the best that can be achieved
for a hub that itself belongs to the Euclidean plane would be distance
$\sqrt{1/2}$.

Now consider a different and simpler optimization criterion: instead of
minimizing distortion, suppose that we want to find a non-contractive mapping
from a given metric space $(X,d)$ to a star network that minimizes the total
or average distance of the points to the hub. That is, for each point $x_{i}$
of the input metric space, we should choose its distance $h_{i}$ to the hub,
such that the mapping from $(X,d)$ to the resulting star metric space is non-
contracting (for every $i$ and $j$, $d(i,j)\leq h_{i}+h_{j}$) and such that we
minimize the sum of the chosen numbers $h_{i}$. Again, this can be viewed as
choosing an optimal point from the tight span of $(X,d)$, but with a different
optimization criterion.

Figure 4: Transforming a metric space $(X,d)$ (upper left) into the modified
metric $(X,d^{*})$ (upper right), finding non-overlapping balls with maximum
sum of radii (lower right), and transforming the radii back into hub distances
(lower left) produces an embedding of $(X,d)$ into a star metric minimizing
the average distance from the hub.

As we now observe, this problem of choosing the star with the minimum average
hub distance can be represented as a problem of minimizing the sum of radii of
non-overlapping disks, for a different metric space $(X,d^{*})$ on the same
set of input points. To see this, let $D$ be the diameter of space $(X,d)$,
the distance between its two farthest points. Note that, in the optimal star
metric space, all hub distances will be at most $D$, for otherwise the
solution would not be minimal. We use $D$ to define the new distance
$d^{*}(x,y)=2D-d(x,y)$. All new distances are between $D$ and $2D$, so this
obeys the triangle inequality and the other requirements of a metric.
Similarly, we convert hub distances $h_{i}$ to disk radii $r_{i}$ by the
formula $r_{i}=D-h_{i}$ or equivalently $h_{i}=D-r_{i}$. Then two hub
distances $h_{i}$ and $h_{j}$ obey the non-contractive inequality $d(i,j)\leq
h_{i}+h_{j}$ if and only if the corresponding two radii $r_{i}$ and $r_{j}$
obey the non-overlapping ball inequality $r_{i}+r_{j}\leq d^{*}(x_{i},x_{j})$.
Therefore, if we find a system of non-overlapping balls in $(X,d^{*})$
minimizing the sum of radii of the balls, this will automatically also give us
a system of hub distances $h_{i}$ defining a non-contractive mapping of
$(X,d)$ to a star network that minimizes the average distance of a point to
the hub. These transformations are illustrated in Figure 4.

If we use these transformations together with the algorithms that we describe
in this paper, we can find the optimal star network embedding for any
$n$-point metric space, minimizing the average distance to the hub, in time
$O(n^{3})$. This compares favorably with the $O(n^{3}\log^{2}n)$ time of
Eppstein and Wortman [10] for finding the star network embedding that
minimizes the distortion.

##  3 Equivalence to cycle cover

Although our radius-sum optimization problem is not equivalent to geometric
matching, we show in this section that it is equivalent by linear programming
duality to a related problem of finding a cycle cover in an associated
geometric graph. We will use this equivalence in formulating algorithms for
both problems.

Let $p_{i}$ ($i=0,\dots n-1$) be a set of points in a metric space, with
distances $d(p_{i},p_{j})$. The problem of finding non-overlapping balls
$(p_{i},r_{i})$ that maximize the sum of the radii can be expressed as a
linear program:

$\text{maximize~{}}\sum r_{i}$

subject to the inequality constraints

$\displaystyle\forall i:\quad$ $\displaystyle r_{i}\geq 0$
$\displaystyle\forall i,j:\quad$ $\displaystyle r_{i}+r_{j}\leq
d(p_{i},p_{j}).$

By linear programming duality [25, Ch. 12] this has the same value as the dual
linear program:

$\text{minimize~{}}\sum w_{ij}d(p_{i},p_{j})$

subject to the inequality constraints

$\displaystyle\forall i,j:\quad$ $\displaystyle w_{ij}\geq 0$
$\displaystyle\forall i:\quad$ $\displaystyle\sum_{j}w_{ij}\geq 1$

That is, we must find non-negative weights $w_{ij}$ for the edges of a
complete geometric graph such that each vertex has incident edge weights
totaling at least one, and minimizing the weighted sum of edge lengths.

######  Lemma .

For distances obeying the triangle inequality, there is an optimal solution to
the dual linear program described above in which the incident edge weights at
each vertex total exactly one.

###### Proof.

Define the excess of a vertex to be the total weight of its incident edges
minus one. If any edge has weight greater than one, its weight can be reduced
to exactly one, improving the quality of the solution without changing its
feasibility. Otherwise, if any vertex has positive excess, we can subtract an
equal positive weight from two of its incident edges and add the same weight
to the edge between the other endpoints of these edges. By the triangle
inequality, this change does not worsen the solution, and again it does not
change its feasibility. By making such changes we can reduce the total excess,
maintaining feasibility, until eventually all vertices have excess zero. ∎

This dual program (either in the form first given above, or with the
restriction that the weights at each vertex sum to exactly one according to
section 3) is the linear programming relaxation of the problem of finding a
minimum weight perfect matching in the complete geometric graph. Such a
relaxation is primarily used for bipartite graphs, for which it is exact [25,
Ex. 12.7]. For non-bipartite graphs such as the complete graph, it has a half-
integral optimal solution in which all weights $w_{ij}$ belong to
$\\{0,1/2,1\\}$ [25, Ex. 14.8]. Doubling the weights to make them integers,
and interpreting each doubled weight as an edge multiplicity, gives us a
combinatorial description of the dual solution as a multiset of edges in which
(by section 3) each vertex has degree two. That is, we have the following
result:

######  Definition .

In a weighted graph $G$, a _cycle cover_ is a multigraph that uses only edges
of $G$ (possibly using some edges twice), and has exactly two edges incident
to every vertex of $G$. The _weight_ of the cycle cover is the sum of weights
of its edges (counted with multiplicity). A _minimum cycle cover_ is a cycle
cover that has the minimum possible weight among all cycle covers of $G$.

######  Theorem 1.

The maximum sum of radii of nonoverlapping balls, centered at points $p_{i}$
of a metric space, equals half of the weight of a minimum cycle cover of the
complete geometric graph on the points $p_{i}$ with edge lengths equal to the
distances between the edge endpoints.

###### Proof.

This follows immediately from the facts that, as with any linear program, the
linear program expressing the problem of maximizing the sum of radii (the
maximization problem at the start of this section) has the same value as its
dual minimization problem, from the known fact that the specific linear
program given by this dual minimization problem has a half-integer solution,
and from section 3 which allows us to interpret this half-integer solution as
a choice of two incident edges at each vertex. ∎

For example, the example of three collinear points, given in subsection 2.1
and Figure 3 as an input that can cause one ball to have zero radius, has a
single three-point cycle as its minimum cycle cover. The length of this cycle,
the sum of the three distances between the points, equals twice the sum of the
radii of the optimal balls (zero for the middle point, and the distance to the
middle point for the two outer points).

Figure 5: Eliminating a long even cycle (here, a rhombus) in the minimum cycle
cover by splitting it into 2-cycles with no greater length.

There always exists a minimum cycle cover in which all cycles of more than two
edges have odd length. To see this, observe that any long even cycle $C$ can
be partitioned into two disjoint matchings, and that at least one of these
matchings gives a covering of the same vertices as $C$ by $2$-cycles whose
total weight is at most that of $C$ (Figure 5). Therefore, from now on we will
assume that our cycle covers have no long even cycles.

Theorem 1 also gives us an easy-to-test optimality condition for the maximum
sum of radii problem, that applies to inputs in general position (without
extra touching balls beyond the ones required by the solution):

######  Definition .

The _touching graph_ of a family of nonoverlapping balls has a vertex for each
ball and an edge connecting each two balls that touch.

######  Corollary .

Suppose that no two balls in a given family of balls overlap, and that each
connected component of the touching graph of the balls is an odd cycle or an
isolated edge. Then this family of balls has the maximum sum of radii of any
family with the same centers.

###### Proof.

The touching graph of the balls (viewing each isolated edge as a $2$-cycle)
gives a cycle cover of length half the sum of radii of the given balls. By
Theorem 1 no cycle cover can be shorter and no system of balls with the same
centers can have a larger sum of radii. ∎

In linear programming terms, this corollary is an instance of _complementary
slackness_ , a more general condition for testing whether a simultaneous
primal and dual solution are both optimal. Figure 6 shows a family of disks
meeting the conditions of the corollary, together with their touching graph.

Figure 6: Nonoverlapping disks that (by section 3) maximize the sum of radii
for their centers.

##  4 Cycle covers from matchings

Theorem 1, giving an equivalence between maximizing the sum of radii and
minimizing the total length of a cycle cover, does not yet give us a
combinatorial algorithm for either of these two equivalent problems. And
because it only gives us an equivalence of the optimum values of the two
problems, it does not help us compute the actual radii in an optimum
collection of disks. Therefore, we still need to show how to solve both
problems combinatorially, in a way that provides for us the solution itself
and not just the solution value.

We begin with the cycle cover. As we show below, it can be transformed into a
more familiar matching problem using the notion of a bipartite double cover.

######  Definition .

The _bipartite double cover_ $2G$ of a graph $G$ is the tensor product
$G\times K_{2}$, a bipartite graph with two copies of each vertex of $G$ (one
of each color) and two copies of each edge of $G$ (one for each pair of
oppositely-colored copies of the endpoint of the edge). If $G$ is weighted, we
use the same edge weights in $G$ and $2G$.

######  Lemma .

Every perfect matching in $2G$ corresponds (under the mapping that takes each
vertex of $2G$ to the corresponding vertex in $G$) to a cycle vertex cover
with equal total length in $G$. Every cycle vertex cover in $G$ comes from a
perfect matching in $2G$ in this way.

###### Proof.

To convert a matching in $2G$ to a cycle cover in $G$, we replace every edge
of the matching by the corresponding edge in $G$. Each vertex in $G$ belongs
to two such edges (one for the matched edge of each of its two copies in $2G$)
so the resulting multiset of edges has degree two at every vertex of $G$.
I.e., it is a cycle cover.

To convert a cycle cover in $G$ to a matching in $2G$ we orient each cycle of
the cycle cover consistently (choosing one of the two possible orientations of
each cycle, arbitrarily). Then, in $2G$, we match each vertex in one copy of
$G$ to its outgoing neighbor in the other copy of $G$.

These two transformations clearly preserve the solution value, proving the
lemma. ∎

See Figure 7 for an example. Each cycle of more than two vertices in $G$ has
two different representations as a set of matched edges in $2G$, but this
ambiguity is not a problem. We can find a minimum cycle vertex cover in $G$ by
finding a minimum weight perfect matching in $2G$.

A minimum weight perfect matching on a bipartite graph of $n$ vertices and $m$
edges can be found in time $O(mn+n^{2}\log n)$ [11]. (Better times are known
when the weights are small integers [7].) By solving this problem on the
doubled complete geometric graph, the optimal sum of radii can be computed in
time $O(n^{3})$. However, this is still slower than we would like for
Euclidean spaces, and does not yet tell us the radii of the individual balls.

Figure 7: The correspondence between a cycle cover of a graph $G$ (left, in
this case, a complete graph $K_{4}$) and a matching in $2G$ (right).

##  5 From cycle covers to balls

Each odd cycle in a minimum cycle cover of the complete geometric graph
corresponds to a unique system of balls that maximizes the sum of radii for
those points. Let the cycle have vertices $p_{0},p_{1},\dots p_{k-1}$ and edge
lengths $\ell_{i}=d(p_{i},p_{i+1\bmod k})$. Then for any $j$ with $0\leq j<i$
set the radius $r_{j}$ of point $p_{j}$ to be

$r_{j}=\sum_{i}\frac{\pm\ell_{i}}{2},$

where we choose the signs in this sum so that the two edges adjacent to
$p_{j}$ have positive sign and every other point $p_{i}$ is incident to two
edges with opposite signs from each other, as depicted in Figure 8.

Figure 8: Converting the edge lengths of an odd cycle to radii of non-
overlapping disks.

######  Lemma .

The radius assignment given above gives the unique set of balls with maximum
sum of radii centered at the vertices of an odd cycle in a minimum cycle
cover.

###### Proof.

For every edge $e$ of the cycle, the two radii of the disks centered at the
endpoints of $e$ have expressions in which the terms for $e$ both have
positive sign. However, these same two expressions differ in their choice of
the signs for every other edge. Therefore, when we sum these two radii, every
term except the ones for the lengths of $e$ cancels, and the sum of the two
radii is exactly $e$ Thus, the two balls centered at the endpoints of $e$
touch but do not overlap, and the sum of all the radii is half the length of
the cycle, as desired. We have not yet shown that balls for non-consecutive
cycle vertices are non-overlapping, but that will follow from uniqueness:
there is no other optimal solution, so this must be a valid solution.

If any other system of non-overlapping balls wth the same centers includes a
consecutive pair that do not touch, then its sum of radii cannot be as large.
For, in any system of nonoverlapping balls, the length of each cycle edge
equals the sum of radii of nonoverlapping balls at its endpoints, plus a non-
negative “gap” equal to the difference between the length and the sum of
radii. So if any gap is positive, the doubled sum of radii will fall short of
the cycle length by that gap.

No other system of non-overlapping balls with the same centers can have each
consecutive pair touching. For, increasing the radius of any one ball
(relative to the solution given above) would cause a chain of alternating
decreases and increases of the radii of all the other balls, of equal
magnitude, leading to an inconsistency because the cycle length is odd and the
changes in radius cannot strictly alternate in sign. ∎

The $2$-cycles in the cycle cover cause us more trouble, because their radii
may be constrained by other nearby points but are not in general uniquely
determined. For instance, in Figure 9, each of the three pairs of touching
circles must have nearly-equal radii to avoid overlaps with nearby circles.

Figure 9: Pairs of disks corresponding to $2$-cycles in the minimum cycle
cover may constrain the radii of other nearby pairs of disks. In this example,
each of the three pairs of touching circles must have nearly-equal radii to
avoid overlaps with nearby circles.

To describe how we transform a minimum weight perfect matching problem on
$2K_{n}$ into a feasible system of balls with maximum sum of radii, we need to
look more deeply into the details of the Hungarian algorithm for minimum
weight perfect matching. The version of this algorithm described by Tarjan
[22, Secs. 8.4 & 9.1] maintains a matching (initially empty) which it uses to
orient the graph. Unmatched edges are directed from one side of the
bipartition to the other, and matched edges are directed in the opposite
direction. In each iteration the algorithm finds a minimum-cost alternating
path between two unmatched vertices; here, the cost of a path is the sum of
lengths of unmatched edges, minus the sum of lengths of matched edges. It uses
this path to extend the matching by one edge.

Because matched edges are subtracted from the path length, each path search
involves negative-weight edges. However, the algorithm also maintains a system
of non-negative weights that are equivalent (in the sense of having the same
shortest paths), allowing Dijkstra’s algorithm to be used, and adjusts these
weights to keep them non-negative after each iteration. Each iteration
increases the number of matched edges, so there are $O(n)$ iterations. The
time per iteration can be bounded by the time for Dijkstra’s algorithm. Using
Fibonacci heaps this gives a total running time of $O(mn+n^{2}\log n)$ [11].

To adjust edge weights we maintain dual variables: a real number for each
vertex of the bipartite graph. We subtract these variables from the length of
each incident unmatched edge, and add these variables to the (negated) length
of each incident matched edge. In order for the adjusted edge weights to be
non-negative, the algorithm maintains an invariant that the length of each
unmatched edge is at least the sum of the dual variables at its endpoints, and
each matched edge length equals the sum of the dual variables at its
endpoints.

Recall that the graph $2K_{n}$ to which we apply this algorithm has two
vertices for each input point $p_{i}$, one of each color (say red and blue).
We may visualize the dual variables at these two vertices as two balls (again,
one red and one blue) both centered at point $p_{i}$. Then the invariants
maintained by the Hungarian algorithm can be expressed in our terms as stating
that pairs of balls of opposite colors cannot overlap unless they have the
same center, and that each matched edge comes from a touching pair of
oppositely colored balls (Figure 10). However, this visualization may be
somewhat misleading, as we allow the dual variables to be negative.

Figure 10: The dual variables for the Hungarian matching algorithm on the
bipartite graph $2K_{n}$ can be visualized as red and blue balls with no
bichromatic overlaps, touching at each matched pair.

######  Lemma .

If we are given the dual variables found by the Hungarian algorithm for
minimum weight perfect matching in the bipartite graph $2K_{n}$, we can
construct from them a system of nonoverlapping balls with maximum sum of
radii, in linear time.

###### Proof.

Let the two dual variables for point $p_{i}$ be $a_{i}$ and $b_{i}$, and
calculate the radius of the ball centered at $p_{i}$ to be the average
$r_{i}=(a_{i}+b_{i})/2$ of the two dual variables. For each $i$ and $j$, we
have (by assumption) $a_{i}+b_{j}\leq d(p_{i},p_{j})$ and $b_{i}+a_{j}\leq
d(p_{i},p_{j})$, and averaging these two inequalities gives us that
$r_{i}+r_{j}\leq d(p_{i},p_{j})$. That is, the balls with radius $r_{i}$ are
non-overlapping.

Each edge of the cycle cover of $K_{n}$ has two balls of opposite color at its
endpoints whose radii sum to the edge length. For each cycle of the cycle
cover of $K_{n}$, each ball centered at a cycle vertex will contribute to this
equality for exactly one of the two incident cycle edges, so the sum of the
radii of the blue and red balls centered at the cycle vertices equals the
length of the cycle. Therefore, when we average these red and blue radii to
give the radii $r_{i}$ of a single system of balls, the sum of the averaged
radii equals half the length of the cycle. Thus, over the whole input, the sum
of all the radii $r_{i}$ equals half the length of the cycle cover. The
optimality of this system of nonoverlapping balls then follows from Theorem 1.
∎

Combining section 4, section 5, and the known algorithms for minimum weight
perfect matching in bipartite graphs gives us the following result.

######  Theorem 2.

Given $n$ points in an arbitrary metric space, we can construct a system of
nonoverlapping balls centered at those points, maximizing the sum of radii of
the balls, in time $O(n^{3})$.

###### Proof.

We construct the bipartite graph $2K_{n}$, weighted by the distances between
points, and find a minimum-weight perfect matching in $2K_{n}$ (or
equivalently by section 4 a cycle vertex cover in $K_{n}$), by using an
algorithm for minimum weight perfect matching that also provides the values of
the dual variables (as the Hungarian algorithm does). We then apply section 5
to convert these dual variable values to an optimal system of radii. ∎

##  6 Euclidean speedup

To speed up the search for balls with maximum sum of radii in low-dimensional
Euclidean spaces, we replace the complete geometric graph of the previous
sections with a sparser graph that contains the optimal cover and behaves the
same with respect to testing whether systems of balls are nonoverlapping.

###  6.1 The nearest-neighbor overlap graph

######  Definition .

Given a system of points $p_{i}$ in a metric space, let the _nearest neighbor
distance_ $\delta_{i}$ of each point $p_{i}$ be

$\delta_{i}=\min_{j\neq i}d(p_{i},p_{j}).$

In Euclidean spaces of bounded dimension, these distances can be found for all
points in total time $O(n\log n)$ [3, 24].

######  Definition .

Define the _nearest-neighbor overlap graph_ $\mathcal{N}$ (Figure 11) to be
the graph whose vertices are the points $p_{i}$, with an edge between two
points $p_{i}$ and $p_{j}$ when the balls with radius equal to the nearest
neighbor distance overlap or touch. That is, $p_{i}$ and $p_{j}$ are adjacent
exactly when

$d(p_{i},p_{j})\leq\delta_{i}+\delta_{j}.$

Figure 11: The nearest neighbor balls of a set of points and their nearest
neighbor overlap graph $\mathcal{N}$.

######  Lemma .

Each edge of the minimum cycle cover of the complete geometric graph belongs
to $\mathcal{N}$.

###### Proof.

Each edge of the cycle cover can be represented by touching balls in a system
of nonoverlapping balls. Because each ball centered at $p_{i}$ in any system
of nonoverlapping balls centered at the given points has radius at most
$\delta_{i}$ (else it would overlap the ball of the nearest neighbor), two
balls of such a system can touch only if the larger balls with radius
$\delta_{i}$ touch or overlap. ∎

######  Lemma .

A system of balls centered at the points $p_{i}$ is nonoverlapping if and only
if, for each edge of $\mathcal{N}$, the two endpoints of the edge are
nonoverlapping.

###### Proof.

A ball centered at $p_{i}$ with radius greater than $\delta_{i}$ overlaps the
ball centered at the nearest neighbor of $p_{i}$, which is necessarily
adjacent to $p_{i}$ in $\mathcal{N}$. Otherwise, if all balls have radius at
most $\delta_{i}$, then they can overlap only if their centers are adjacent in
$\mathcal{N}$. ∎

By subsection 6.1 we can find our minimum cycle cover by applying a weighted
matching algorithm to $2\mathcal{N}$. The system of radii formed by averaging
the dual variables of the matching in $2\mathcal{N}$ will give us an optimal
set of nonoverlapping balls, just as it did in $2K_{n}$. For, the averaged
radii will achieve the optimal value and will have no overlaps that are
detected by the edges of $\mathcal{N}$; by subsection 6.1 there can be no
other overlaps. It remains to show that we can find a matching (and its dual
variables) quickly in this graph. For this we will need some geometry, since
in arbitrary metric spaces $\mathcal{N}$ can be complete.

###  6.2 Weighted matching in separated bipartite graphs

In this section we prove the existence of an efficient minimum weight perfect
matching algorithm for bipartite graphs obeying a separator theorem. We will
later show that the nearest-neighbor graph $\mathcal{N}$ and its bipartite
double cover obey such a theorem in low-dimensional Euclidean spaces, allowing
us to apply this algorithm.

######  Definition .

Let $s$ be a real number in the range $0<s<1$. Define an $s$-_separator_ of an
$n$-vertex graph to be a subset $S$ of the vertices the removal of which
allows the graph to be partitioned into two subgraphs, disconnected from each
other and each having at most $sn$ vertices. When $s$ is a constant (depending
on some general family of graphs but not the specific graph in that family) we
will omit it and call $S$ a _separator_.

An efficient algorithm for weighted matching in graphs with small separators
was already given by Lipton and Tarjan [19, Sec. 7]. They directly consider
only planar graphs, but their method would also work for other graphs that
have separator hierarchies. However, they consider a different variant of
matching, maximum weight matching on non-bipartite graphs, and are not
explicit in how they maintain and update the analogue for that problem of the
dual variables that we need. They also do not take advantage of more recent
advances in shortest path algorithms for graphs with separators [13]. Here we
describe an algorithm based on similar ideas to the Lipton–Tarjan algorithm
that maintains the dual variables that we need and is faster by a logarithmic
factor.

Suppose we are given a weighted graph $G$ with a separator $S$ (part of a
separator hierarchy for $G$), such that $G\setminus S$ can be partitioned into
two subgraphs $L$ and $R$, disconnected from each other, with $L$ and $R$ both
having a number of vertices that is smaller by a constant factor. For such a
graph, Lipton and Tarjan [19, Sec. 7] describe a divide-and-conquer algorithm
for maximum-weight matching (an equivalent problem to minimum-weight perfect
matching) that recursively matches $L$ and $R$, and then adds the separator
vertices one at a time to the union of $L$ and $R$. When each separator vertex
is added, the maximum matching may change along a single alternating path.

Lipton and Tarjan write that “given a suitable representaton”, each
alternating path can be found in the time for a single application of
Dijkstra’s algorithm. This, in turn, allows the most expensive term in the
running time of their algorithm to be bounded by the product of the time for
finding each path with the number of separator vertices. The “suitable
representation” that they refer to is made more complicated by the possibility
that $G$ might not be bipartite, but in the bipartite case it is essentially
the system of dual variables that we need. However, their handling of these
variables is lacking in detail, so we provide here a more explicit divide and
conquer matching algorithm, with the following differences:

  * •

We formulate the problem as finding a minimum-weight maximum-cardinality
matching in a graph with non-negative edge weights. This slight generalization
of minimum-weight perfect matching allows us to perform recursive calls
without needing to guarantee that each subproblem has a perfect matching, and
simplifies the definition of the dual variables relative to their definition
for maximum-weight matching.

  * •

We recurse on the induced subgraphs $G[S\cup L]$ and $G[S\cup R]$ rather than
on $L$ and $R$. This inclusion of the separator vertices in the recursive
subproblems results in the calculation of two inconsistent systems of dual
variables for these vertices, which must be reconciled when the subproblem
solutions are combined.

  * •

We compute shortest paths using an algorithm of Henzinger et al. [13] rather
than by Dijkstra’s algorithm. This algorithm takes linear time per path, after
a separator hierarchy has already been computed, saving a logarithmic time
factor compared to Dijkstra.

  * •

We consider only bipartite graphs rather than the more general case of
weighted matching in arbitrary graphs.

Because we recurse on $G[S\cup L]$ and $G[S\cup R]$, we us a separator
hierarchy in which these two subgraphs (rather than $L$ and $R$) are the ones
that are recursively subdivided. In more detail, our algorithm performs the
following steps:

  1. 1.

Recursively compute a minimum-weight maximum-cardinality matching, and its
corresponding system of dual variables, for the two subgraphs $G[S\cup L]$ and
$G[S\cup R]$. Recall that the dual variables are numbers associated with each
vertex such that, for each unmatched edge, the two numbers at its endpoints
sum to at most the edge length, and for each matched edge they sum to exactly
the edge length.

  2. 2.

Combine the two subproblems to give a single matching (with fewer than the
maximum number of matched edges) and valid system of dual variables for the
whole graph $G$. To do so, for each vertex $s$ in $S$, we choose the dual
variable for $s$ in $G$ to be the minimum of the two values computed for $s$
in $G[L\cup S]$ and in $G[R\cup S]$, and we remove the matched edge incident
to this vertex in one of the two subproblems, the subproblem that did not
supply the minimum dual variable value. (If the two subproblems have equal
dual variables at $s$, we choose arbitrarily which of the two matched edges to
remove.)

  3. 3.

Repeat the following steps until no more alternating paths can be found:

     * •

Modify the original weight of each edge by subtracting the dual variables at
its endpoints.

     * •

Add an artificial source vertex adjacent by zero-length edges to all unmatched
vertices on one side of the bipartition. Orient the unmatched edges of the
graph from the source side to the opposite side, and the matched edges in the
other direction.

     * •

Use the algorithm of Henzinger et al. starting from this source vertex to find
the alternating path of minimum modified weight between two unmatched
vertices.

     * •

Augment the matching using the alternating path, trading matched and unmatched
edges along the path.

     * •

Use the distances from the source vertex found by the path search to update
the dual variables so that they remain valid.

To update the dual variables, decrease each dual variable on the source side
of the bipartition by its distance from the source, and increase each dual
variable on the opposite side by its distance from the source. The sum of dual
variables at the endpoints of each matched edge remains unchanged, the sum of
variables on any shortest path edge increases to equal the edge length, and
the sum of variables on any other edge remains at most equal to the edge
length (else that edge would have supplied a shorter path).

The result of this algorithm is a maximum-cardinality matching on $G$, and a
feasible system of dual variables, whose existence ensures that the matching
has minimum weight (as any alternating cycle would have non-negative modified
weight and therefore non-negative total cost). This proves the following
result:

######  Lemma .

Let $G$ be a bipartite graph given together with a separator hierarchy in
which each subgraph of $p$ vertices has a separator of size $O(p^{c})$, for
some constant $c$ with $0<c<1$. Then one can find both a minimum weight
perfect matching of $G$, and the dual variables of the matching, in time
$O(n^{1+c})$.

###### Proof.

The union of the matchings in $G[S\cup L]$ and $G[S\cup R]$ is at least as
large as the single maximum matching in $G$, so the reconciliation process (in
which we remove at most $|S|$ edges from the union of the matchings) produces
a single matching that differs in cardinality from the maximum by at most
$|S|$. Therefore, the inner loop of the algorithm involves at most $|S|$
linear-time path searches. The time analysis of subsection 6.2 follows
straightforwardly as a divide-and-conquer recurrence. ∎

###  6.3 Low-ply neighborhood systems

To show that the nearest-neighborhood graph has small separators, and hence
that we can quickly solve weighted matching problems on it, we need some
definitions from Teng [23].

######  Definition .

A _$k$ -neighborhood system_ is a system of balls such that each ball contains
at most $k$ centers of balls in the system.

Clearly, the nearest-neighbor balls from which $\mathcal{N}$ is defined form a
$1$-neighborhood system, because each contains only its own center.

######  Definition .

The _ply_ of a system of balls is the maximum number of balls that have a
common point of intersection. The _kissing number_ $\tau_{d}$ of
$d$-dimensional Euclidean space is the maximum number of unit-radius balls
that can touch a single central unit ball without overlapping each other.

The kissing number is known to depend singly-exponentially on the dimension,
but for us it is more important that it is constant for any constant
dimension.

######  Lemma  (Teng [23], Lemma 3.2).

Let $k$ and $d$ be fixed constants. Then every $k$-neighborhood system in
$d$-dimensional Euclidean space has ply at most $\tau_{d}k=O(1)$.

The following definition is not from Teng (who uses different terminology) but
is standard in graph theory.

######  Definition .

The _degeneracy_ of a graph $G$ is the smallest value $d$ such that every
subgraph of $G$ contains a vertex of degree at most $d$.

######  Lemma  (Teng [23], Theorem 4.2).

For any $k$-neighborhood system in a Euclidean space of bounded dimension, the
intersection graph of the balls has bounded degeneracy.

In particular, this implies that $\mathcal{N}$ has $O(n)$ edges. This
immediately gives a speedup from $O(n^{3})$ to $O(n^{2}\log n)$, by applying
the $O(mn+n^{2}\log n)$ time bound for matching on the sparse bipartite graph
$2\mathcal{N}$ and using the degeneracy bound to substitute $m\leq
n\Delta=O(n)$. But we can do better using more graph structure, a separator
theorem for overlap graphs.

######  Lemma  (Teng [23], Theorem 5.3).

For any system of balls of bounded ply in a Euclidean space of bounded
dimension $d$, the intersection graph of the balls has separators of size
$O(n^{1-1/d})$

The following deterministic algorithmic version of this separator theorem is
from [9]:

######  Lemma  (Eppstein et al. [9], Theorems 5.1 and 5.8).

For any system of balls of bounded ply in a Euclidean space of bounded
dimension $d$, the intersection graph of the balls has separators of size
$O(n^{1-1/d})$ that can be found deterministically in linear time from the
balls. By recursively constructing a hierarchy of separators, the overlap
graph can be constructed from the balls in time $O(n\log n)$.

Another way of expressing the existence of small separators (of size a
fractional power of $n$) is that overlap graphs are graphs of _polynomial
expansion_ [8, 20]. We remark that any hierarchy of separators for
$\mathcal{N}$ can be transformed into a hierarchy of separators for
$2\mathcal{N}$, simply by including both copies of each vertex of a separator
for $\mathcal{N}$ in the corresponding separator for $2\mathcal{N}$.

###  6.4 Convex distance functions

A _convex distance function_ may be defined from any centrally symmetric
convex body $B$ on $d$-dimensional Euclidean space by defining the distance
between any two points $p$ and $q$ to be the infimum of the positive numbers
$s$ such that a copy of $B$, centered at $p$ and expanded by a factor of $s$,
contains $q$. The distance defined in this way obeys the axioms of a metric
space, and gives the Euclidean space the structure of a normed vector space.
Conversely, the distances in any bounded dimensional normed vector space can
be defined in this way from the unit ball of the norm. Much, but not all, of
the same theory for Euclidean balls goes through in these spaces. In
particular, we have:

######  Lemma  (Teng [23], Lemma 7.2).

Let $k$ be a fixed constant, fix a convex distance function on a normed metric
space of fixed dimension $d$, and let $\tau$ be the kissing number for that
space (the maximum number of non-overlapping unit balls that can
simultaneously touch a single central unit ball). Then Every $k$-neighborhood
system in the space has ply at most $\tau_{d}k\leq(3^{d}-1)k=O(1)$.

######  Lemma  (Teng [23], Lemma 7.4).

Let $k$ be a fixed constant, and fix a convex distance function on a normed
metric space of fixed dimension $d$. Then the intersection graph of any
$k$-neighborhood system has bounded degeneracy.

######  Lemma  (Teng [23], Theorem 7.1).

Let $k$ be a fixed constant, and fix a convex distance function on a normed
metric space of fixed dimension $d$. Then the intersection graph of any
$k$-neighborhood system has separators of size $O(n^{1-1/d})$ that can be
found by a randomized algorithm in linear expected time from the balls. By
recursively constructing a hierarchy of separators, the overlap graph can be
constructed from the balls in randomized expected time $O(n\log n)$.

###  6.5 Optimal disks for low-dimensional Euclidean spaces

Putting the results of this section together, we have our main theorem:

######  Theorem 3.

Given $n$ points in a Euclidean space of constant dimension $d\geq 2$, we can
construct a system of balls centered at the given points, with maximum sum of
radii, in time $O(n^{2-1/d})$. For points in a metric space of dimension $d$
defined by a convex distance function, we can solve the same problem in
randomized expected time $O(n^{2-1/d})$.

###### Proof.

We first find the nearest neighbors of all points. For this, the deterministic
algorithm of Vaidya [24] takes time $O(n\log n)$. It is described by Vaidya
for $L_{p}$ metrics, but depends only on the property that, if the space is
partitioned into a grid of congruent cubes, then only a bounded number of
cubes can intersect a ball whose radius is the diameter of a single cube. This
is true for any convex distance function, so Vaidya’s algorithm can be used to
find all nearest neighbors for any convex distance function.

We then construct graph $\mathcal{N}$ and its separator hierarchy by
subsection 6.3 (for Euclidean distances) or subsection 6.4 (for convex
distance functions). We construct the bipartite double cover $2\mathcal{N}$
and a system of separators for it that are formed by doubling the separators
in $\mathcal{N}$. Using these separators, we apply the matching algorithm of
subsection 6.2 to $2\mathcal{N}$, giving two dual variables per input point.
Finally, we take the average of these two dual variables to give a single
radius for each point. ∎

##  7 Lower-bounding the radii

As discussed in subsection 2.1, in the map-labeling application of our problem
it is helpful to constrain the radii of the disks or balls to all be at least
$\delta$, for some given value $\delta>0$ (small enough that all
$\delta$-balls are non-overlapping), so that we avoid zero-radius or overly-
small balls in our solution.

This constrained problem, for any given value of $\delta$, can be transformed
to an instance of our unconstrained radius-sum optimization problem by
modifying the input distances, as follows.

  * •

Form a collection of balls of radius $\delta$ around each given point $p_{i}$,
and let $d^{\circ}(p_{i},p_{j})=d(p_{i},p_{j})-2\delta$ be the distance
between each pair of balls. Then $d^{\circ}$ is non-negative and symmetric,
but it might not obey the triangle inequality, which we will need in our
algorithms (in particular, in section 3).

  * •

Let $d^{*}$ be the shortest path distance in a complete geometric graph whose
edge weights are $d^{\circ}$. The distance $d^{*}$ defined in this way is non-
negative, symmetric, and obeys the triangle inequality.

  * •

Find a system of non-negative radii $r_{i}^{*}$ of disjoint metric balls for
$d^{*}$, maximizing the sum of radii.

  * •

Set $r_{i}=r_{i}^{*}+\delta$.

######  Theorem 4.

For $r_{i}$ computed by the steps outlined above, the balls of radius $r_{i}$
around each point $p_{i}$ of the original point set all have radius at least
$\delta$, and have the maximum total radius of any system of balls of radius
at least $\delta$.

###### Proof.

Let $r_{i}^{\circ}$ be a system of radii, obeying the constraints that for
each pair of points $p_{i}$ and $p_{j}$, $r_{i}^{\circ}+r_{j}^{\circ}\leq
d^{\circ}(p_{i},p_{j})$ and maximizing the sum of radii. Then the constraint
that $r_{i}^{\circ}+r_{j}^{\circ}\leq d^{\circ}(p_{i},p_{j})$ is equivalent
(by simple algebraic manipulation) to the constraint that

$(r_{i}^{\circ}+\delta)+(r_{j}^{\circ}+\delta)\geq d(p_{i},p_{j}),$

or in words that the two ball radii derived from $r_{i}^{\circ}$ and
$r_{j}^{\circ}$ form non-overlapping balls around the original points.
Clearly, in addition, these radii are at least $\delta$ and are otherwise
unconstrained. So if we could find the system of radii $r_{i}^{\circ}$, the
radii $r_{i}^{\circ}+\delta$ derived from them would solve the problem stated
in the theorem. It remains to show that each radius $r_{i}^{*}$ found by the
procedure outlined above actually equals $r_{i}^{\circ}$.

To see this, consider the transformation of distances made in going from
$d^{\circ}$ to $d^{*}$. This transformation can only decrease any pairwise
distance; that is, for every $i$ and $j$,

$d^{*}(p_{i},p_{j})\leq d^{\circ}(p_{i},p_{j}).$

For, the distance $d^{\circ}(p_{i},p_{j})$ is the length of the edge
$p_{i}p_{j}$ in the complete weighted graph used to compute $d^{*}$, and this
edge forms one of the paths in the shortest path distance computation. Because
$d^{*}$ has distances that are no greater than $d^{\circ}$, the radii $r^{*}$
must also be no greater than the corresponding values of $r^{\circ}$.

To complete the proof, we show that the system of radii $r^{\circ}$
automatically obey the constraints $r_{i}^{\circ}+r_{j}^{\circ}\leq
d^{*}(p_{i},p_{j})$ used to determine the radii $r^{*}$. Since both
$r^{\circ}$ and $r^{*}$ are defined by optimization problems with the same
objective function but with tighter constraints on $r^{*}$, and since these
tighter constraints are already obeyed by $r^{\circ}$, it will follow that the
optimization problems defining $r^{\circ}$ and $r^{*}$ have identical solution
sets. Therefore, consider any two points $p_{i}$ and $p_{j}$, and their radii
$r_{i}^{\circ}$ and $r_{j}^{\circ}$; we must show that
$r_{i}^{\circ}+r_{j}^{\circ}\leq d^{*}(p_{i},p_{j})$. To see this, consider
the shortest path from $p_{i}$ to $p_{j}$ in the complete weighted graph from
which we defined $d^{*}$, and let $S$ be the sum, over all the edges in this
path, of the radii $r^{\circ}$ at the two endpoints of this edge. Then each
edge’s contribution to the sum is less than its length, for otherwise the two
circles of radius $r^{\circ}$ at its endpoints would violate the constraint
that their sum is at most $d^{\circ}$. Therefore, $S$ is less than the total
length of the path, which is $d^{*}(p_{i},p_{j})$. But both $r_{i}^{\circ}$
and $r_{j}^{\circ}$ contribute to $S$ in the terms for the first and last edge
to the path, and all the other contributions are non-negative. Therefore,
$r_{i}^{\circ}+r_{j}^{\circ}\leq S\leq d^{*}(pi,p_{j})$, as desired. ∎

We leave open for future research the problem of combining this technique for
lower-bounding the radii of systems of balls with our speedups for Euclidean
radius-sum optimization.

## References

  * [1] H. Alt, E. M. Arkin, H. Brönnimann, J. Erickson, S. P. Fekete, C. Knauer, J. Lenchner, J. S. B. Mitchell, and K. Whittlesey.  Minimum-cost coverage of point sets by disks.  Proc. 22nd Symp. Computational Geometry (SoCG 2006), pp. 449–458, 2006, doi:10.1145/1137856.1137922, MR2389355.
  * [2] E. W. Chambers, S. P. Fekete, H.-F. Hoffmann, D. Marinakis, J. S. B. Mitchell, V. Srinivasan, U. Stege, and S. Whitesides.  Connecting a set of circles with minimum sum of radii.  Proc. 12th Int. Symp. Algorithms and Data Structures (WADS 2011), pp. 183–194. Springer, LNCS 6844, 2011, doi:10.1007/978-3-642-22300-6_16, MR2863135.
  * [3] K. L. Clarkson.  Fast algorithms for the all nearest neighbors problem.  Proc. 24th IEEE Symp. Foundations of Computer Science (FOCS ’83), pp. 226–232, 1983, doi:10.1109/SFCS.1983.16.
  * [4] E. Cohen and N. Megiddo.  Improved algorithms for linear inequalities with two variables per inequality.  SIAM J. Comput. 23(6):1313–1347, 1994, doi:10.1137/S0097539791256325, MR1303338.
  * [5] W. J. Cook, W. H. Cunningham, W. R. Pulleyblank, and A. Schrijver.  Combinatorial optimization.  Wiley-Interscience Series in Discrete Mathematics and Optimization. John Wiley & Sons, 1998, MR1490579.
  * [6] S. Doddi, M. V. Marathe, A. Mirzaian, B. M. E. Moret, and B. Zhu.  Map labeling and its generalizations.  Proc. 8th ACM-SIAM Symp. on Discrete Algorithms (SODA ’97), pp. 148–157, 1997, MR1447660.
  * [7] R. Duan and H.-H. Su.  A scaling algorithm for maximum weight matching in bipartite graphs.  Proc. 23rd ACM-SIAM Symp. on Discrete Algorithms (SODA ’12), pp. 1413–1424, 2012, MR3205301.
  * [8] Z. Dvořák and S. Norin.  Strongly sublinear separators and polynomial expansion.  Electronic preprint arxiv:1504.04821, 2015.
  * [9] D. Eppstein, G. L. Miller, and S.-H. Teng.  A deterministic linear time algorithm for geometric separators and its applications.  Fundamenta Informaticae 22(4):309–331, 1995, MR1360950.
  * [10] D. Eppstein and K. A. Wortman.  Optimal embedding into star metrics.  Proc. Algorithms and Data Structures Symposium (WADS 2009), pp. 290–301. Springer-Verlag, Lecture Notes in Computer Science 5664, 2009, doi:10.1007/978-3-642-03367-4_26, arXiv:0905.0283, MR2550615.  Winner, best paper award.
  * [11] M. L. Fredman and R. E. Tarjan.  Fibonacci heaps and their uses in improved network optimization algorithms.  J. ACM 34(3):596–615, 1987, doi:10.1145/28869.28874, MR904195.
  * [12] M. Gibson, G. Kanade, E. Krohn, I. A. Pirwani, and K. Varadarajan.  On metric clustering to minimize the sum of radii.  Algorithmica 57(3):484–498, 2010, doi:10.1007/s00453-009-9282-7, MR2609050.
  * [13] M. R. Henzinger, P. Klein, S. Rao, and S. Subramanian.  Faster shortest-path algorithms for planar graphs.  J. Comput. System Sci. 55(1):3–23, 1997, doi:10.1006/jcss.1997.1493, MR1473046.
  * [14] D. S. Hochbaum and J. Naor.  Simple and fast algorithms for linear and integer programs with two variables per inequality.  SIAM J. Comput. 23(6):1179–1192, 1994, doi:10.1137/S0097539793251876, MR1303329.
  * [15] P. Indyk and J. Matoušek.  Low-distortion embeddings of finite metric spaces.  Handbook of Discrete and Computational Geometry, 2nd edition, pp. 177–196. CRC, Discrete Mathematics and its Applications (Boca Raton), 2004, doi:10.1201/9781420035315.ch8.
  * [16] J. R. Isbell.  Six theorems about injective metric spaces.  Comment. Math. Helv. 39:65–76, 1964, doi:10.1007/BF02566944.
  * [17] M. Jiang, S. Bereg, Z. Qin, and B. Zhu.  New bounds on map labeling with circular labels.  Proc. 15th Int. Symp. Algorithms and Computation (ISAAC 2004), pp. 606–617. Springer, LNCS 3341, 2004, doi:10.1007/978-3-540-30551-4_53, MR2158365.
  * [18] M. Jünger and W. Pulleyblank.  New primal and dual matching heuristics.  Algorithmica 13(4):357–380, 1995, doi:10.1007/BF01293485, MR1318310.
  * [19] R. J. Lipton and R. E. Tarjan.  Applications of a planar separator theorem.  SIAM J. Comput. 9(3), 1980, doi:10.1137/0209046, MR584516.
  * [20] J. Nešetřil and P. Ossona de Mendez.  Sparsity: Graphs, Structures, and Algorithms.  Algorithms and Combinatorics 28. Springer, 2012, doi:10.1007/978-3-642-27875-4, MR2920058.
  * [21] T. Strijk and A. Wolff.  Labeling points with circles.  Int. J. Comput. Geom. Appl. 11(02):181–195, 2001, doi:10.1142/s0218195901000444, MR1831031.
  * [22] R. E. Tarjan.  Data Structures and Network Algorithms.  CBMS-NSF Regional Conference Series in Applied Mathematics 44. SIAM, 1983, doi:10.1137/1.9781611970265, MR826534.
  * [23] S.-H. Teng.  Points, spheres, and separators: a unified approach to graph partitioning.  Ph.D. thesis, Carnegie-Mellon Univ., School of Computer Science, 1991, MR2687483.
  * [24] P. M. Vaidya.  An $O(n\log n)$ algorithm for the all-nearest-neighbors problem.  Discrete Comput. Geom. 4(2):101–115, 1989, doi:10.1007/BF02187718, MR973540.
  * [25] V. V. Vazirani.  Approximation Algorithms.  Springer, 2001, doi:10.1007/978-3-662-04565-7, MR1851303.
  * [26] F. Wagner and A. Wolff.  A practical map labeling algorithm.  Comp. Geom. Th. & Appl. 7(5-6):387–404, 1997, doi:10.1016/s0925-7721(96)00007-7, MR1447248.