# Atomic electron shell excitations
in double-$\beta$ decay

M.I. Krivoruchenko†, K.S. Tyrin†, F.F. Karpeshin∗ †National Research Centre
“Kurchatov Institute”
Ploschad’ Akademika Kurchatova 1, 123182 Moscow, Russia
∗D. I. Mendeleev Institute for Metrology (VNIIM)
Moskovsky Ave. 19, 190005 Saint Petersburg, Russia

###### Abstract

The problem of the transition of electron shells of atoms to excited states in
the process of neutrinoless double-$\beta$ decay is investigated. This subject
is crucial for modeling the energy spectrum of $\beta$-electrons, which is
sensitive to the mass and Majorana nature of neutrinos. The dependence of the
obtained results on the atomic number indicates an important role of the
Feinberg–Migdal effect in the electron shell excitations. We report the
overlap amplitudes of the electron shells of the parent atom and the daughter
ion for eleven atoms, the two-neutrino double-$\beta$ decay of which was
observed experimentally. In around one-fourth of the cases where the structure
of the electron shells is inherited from the parent atom, there is a
transition to the ground state or the excited state with the lowest energy.
The de-excitation of the daughter ion in the latter scenario is accompanied by
the emission of photons in the ultraviolet range, which can serve as an
auxiliary signature of double-$\beta$ decay. The average excitation energy of
the electron shells ranges between 300 and 800 eV, with the variance ranging
from $(1.7~{}\mathrm{keV})^{2}$ in calcium to $(14~{}\mathrm{keV})^{2}$ in
uranium.

††journal: JETP Letters

Neutrinoless double-$\beta$ decay (0$\nu$2$\beta$) does not preserve the total
number of leptons and is particularly interesting in the search for departures
from the Standard Model (SM). Similar significance could be found in the quark
sector of SM for processes which violate baryon number conservation, such as
proton decay and neutron-antineutron oscillations [1]. Beyond the SM, any
mechanism of 0$\nu$2$\beta$ decay implies the existence of a Majorana neutrino
mass [2, 3]. In the effective theory, the Majorana neutrino mass, $m_{\nu}$,
is generated by the Weinberg operator of dimension $d=5$ [4]. In the absence
of operators of dimension $d>5$ and symmetry between left and right elementary
fermions, the amplitude of 0$\nu$2$\beta$ decay with light Majorana neutrino
is proportional to $m_{\nu}$.

Experimental searches for 0$\nu$2$\beta$ decay have been actively preformed
for a number of decades. The GERDA collaboration recently obtained a
constraint $m_{\nu}<0.079-0.18$ eV at the confidence level CL = 90% using the
isotope 76Ge [5]. Similar results were obtained by the EXO collaboration [6]
using xenon-136. Restriction on the neutrino Majorana mass $m_{\nu}<0.3-0.9$
eV was also obtained by the collaboration NEMO-3 using molybdenum-100 [7]. The
SuperNEMO experiment is under preparation [8]. CUORE’s experiments using the
isotope 130Te [9, 10] and KamLAND-Zen with liquid xenon-136 [11] are in the
active phase.

The uncertainty of the upper limit on the neutrino mass is due to the accuracy
of calculations of the nuclear part of the process [12, 13, 14].

Experiments to search for 0$\nu$2$\beta$ decay analyse the energy spectrum at
the boundary of the phase space of $\beta$-electrons in order to find a
deviation from the energy spectrum of a more probable two-neutrino
double-$\beta$ decay (2$\nu$2$\beta$). Experimentalists inevitably encounter a
problem that has become widely known in connection with attempts to measure
neutrino mass in tritium beta decay: the daughter atom with a high probability
passes into an excited state. This may be the excited state of a molecule
composed of active target atoms. The atoms themselves experience excitation
due to shake-up and shake-off effects or internal scattering of
$\beta$-electrons. The theory of these processes is developed by Feinberg [39]
and Migdal [40]. Influence of these processes on the spectrum of
$\beta$-electrons are especially noticeable near the spectrum boundary. The
effect increases significantly due to the fact that the spread of residual
excitation energies is almost an order of magnitude higher than the average
value [15]. A similar effect can be expected from the chemical shift [16].

The implications of atom ionization and excitation, first studied in the
context of nuclear physics, are observed in molecular, solid-state systems and
are crucial to the experiments LUX [17], XENON1T [18], and DarkSide-50 [19],
which are designed to detect dark matter particles.

In double-$\beta$ decays, the daughter ion with a high probability occurs in
an excited state [20, 21, 22, 23], which reduces the energy carried away by
$\beta$-electrons. The energy spectrum of $\beta$-electrons in 0$\nu$2$\beta$
decay is a delta function, distorted by atomic effects. This peak is
considered as the 0$\nu$2$\beta$ decay’s signature. The decay realizes a
scenario in which channels with valence electron excitations prevail in
probability, although the average excitation energy, $\mathcal{M}$, and its
variance, $\mathcal{D}$, are essentially saturated by rare electron
excitations from inner atomic orbitals.

In this paper, we estimate the deviations of the $\beta$-electrons energy from
the decay energy, $Q^{*}$, of the 0$\nu$2$\beta$ decay for 11 atoms for which
2$\nu$2$\beta$ decay was experimentally observed.

The binding energy of electrons on the K shell differs from the binding energy
of valence electrons by around three orders of magnitude ($\sim Z^{2}$) in
medium-heavy and heavy atoms, making it difficult to estimate the magnitude of
$\mathcal{M}$ and $\mathcal{D}$ qualitatively. Excitation of valence electrons
with low binding energy obviously dominates the decay probability. However,
the calculations result in unusually large values of average excitation energy
and its variance. Given that the accuracy of calculations in multi-particle
problems is also limited, the paper considers several approaches, including
the Thomas-Fermi (TF) model [24], the Thomas-Fermi-Dirac-Weizsacker model
(TFDW) [25, 26, 27, 28, 29], non-relativistic Rutaan-Hartree-Fock (RHF)
formalism [30] and relativistic Dirac-Hartree- Fock formalism (DHF) [31, 32,
33, 34, 35]. When the outcomes are compared, the magnitude of uncertainty in
the parameters of interest can be evaluated.

Each of these approaches has its advantages and limitations. Unlike the TF
model, in the TFDW model the electron density is finite at the nucleus, which
makes it possible to determine the variance within the model. In the RHF
method, the wave functions of orbitals are parametrized analytically, which
makes it possible to find exchange contributions to the variance and other
observables, but the applicability of the method is restricted to light and
medium-heavy atoms. Within the framework of DHF, the basic properties of
atomic electron shells are tabulated in [31, 32, 33] and implemented in the
form of software packages such as Grasp-2018 [34, 35] and RAINE [36, 37].

In what follows, the system of atomic units $\hbar=m=e=1$, $c=137$ is used,
where $m$ is the electron mass, $e$ is the proton charge, $c$ is the speed of
light. Let $\hat{H}_{Z,N}$ be the Hamiltonian of $N$ electrons of an ion with
a nucleus charge $Z$. We denote $|Z,N\rangle$ the ground state and $E_{Z,N}$
the binding energy of the electrons, so that
$\hat{H}_{Z,N}|Z,N\rangle=E_{Z,N}|Z,N\rangle$.

Table 1:  The average excitation energy of the electron shells of the daughter ion and the variance for eleven atoms, the $2\nu 2\beta$ decay of which was observed experimentally. The second column contains the values of the mass difference, $Q$, of the neutral atoms involved in the decay. The fourth column shows the overlap integrals of the electron shells of the parent atom and the twice ionized daughter atom, calculated with the use of the software package Grasp-2018 [34, 35]. The average energy of the electron shell excitations of the daughter ion is shown in the TF, TFDW and DHF models, the upper bound of the variance $\bar{\mathcal{D}}$ is shown in the TF, TFDW models, and the variance $\mathcal{D}$ – in the DHF and RHF models. The values of $\mathcal{M}_{\mathrm{DHF}}$ and $\mathcal{D}_{\mathrm{DHF/a}}$ excluding exchange contributions are obtained using the results of [33] and [32], respectively. $\mathcal{D}_{\mathrm{DHF/b}}$ includes exchange contributions. To calculate $\mathcal{D}_{\mathrm{RHF/a}}$ without and $\mathcal{D}_{\mathrm{RHF/b}}$ with exchange contributions, the wave functions of orbitals in the RHF method are used [30]. The double ionization energy $I_{2}$ [38] is rounded to three significant digits. The predictions of the non-relativistic models TF, TFDW and RHF are limited by the values of the nuclear charge $Z\leq 54$.  Process | $\begin{array}[]{c}Q\\\ \textup{[keV]}\end{array}$ | Ref. | $K_{Z}$ | $\begin{array}[]{c}\mathcal{M}_{\mathrm{TF}}\\\ \textup{[eV]}\end{array}$ | $\begin{array}[]{c}\mathcal{M}_{\mathrm{TFDW}}\\\ \textup{[eV]}\end{array}$ | $\begin{array}[]{c}\mathcal{M}_{\mathrm{DHF}}\\\ \textup{[eV]}\end{array}$ | $\begin{array}[]{c}I_{2}\\\ \textup{[eV]}\end{array}$ | $\begin{array}[]{c}\bar{\mathcal{D}}_{\mathrm{TF}}^{1/2}\\\ \textup{[keV]}\end{array}$ | $\begin{array}[]{c}\bar{\mathcal{D}}_{\mathrm{TFDW}}^{1/2}\\\ \textup{[keV]}\end{array}$ | $\begin{array}[]{c}\mathcal{D}_{\mathrm{DHF/a}}^{1/2}\\\ \textup{[keV]}\end{array}$ | $\begin{array}[]{c}\mathcal{D}_{\mathrm{DHF/b}}^{1/2}\\\ \textup{[keV]}\end{array}$ | $\begin{array}[]{c}\mathcal{D}_{\mathrm{RHF/a}}^{1/2}\\\ \textup{[keV]}\end{array}$ | $\begin{array}[]{c}\mathcal{D}_{\mathrm{RHF/b}}^{1/2}\\\ \textup{[keV]}\end{array}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
${}_{\;\,20}^{\;\,48}\mathrm{Ca}$ | $\rightarrow$ | ${}_{\;\,22}^{\;\,48}\mathrm{Ti}$ | 4267.98(32) | [42] | 0.466 | 335 | 247 | 299 | 20.4 | 1.25 | 2.43 | 1.70 | 1.65 | 1.66 | 1.61
${}_{\;\,32}^{\;\,76}\mathrm{Ge}$ | $\rightarrow$ | ${}_{\;\,34}^{\;\,76}\mathrm{Se}$ | 2039.006(50) | [43] | 0.575 | 383 | 246 | 369 | 30.9 | 2.16 | 3.92 | 2.88 | 2.77 | 2.72 | 2.62
|  |  | 2039.061(7) | [5] |  |  |  |  |  |  |  |  |  |  |
${}_{\;\,34}^{\;\,82}\mathrm{Se}\,$ | $\rightarrow$ | ${}_{\;\,36}^{\;\,82}\mathrm{Kr}$ | 2997.9(3) | [44] | 0.597 | 384 | 238 | 377 | 38.4 | 2.31 | 4.17 | 3.09 | 2.97 | 2.90 | 2.79
${}_{\;\,40}^{\;\,96}\mathrm{Zr}\,$ | $\rightarrow$ | ${}_{\;\,42}^{\;\,96}\mathrm{Mo}$ | 3356.097(86) | [45] | 0.518 | 422 | 246 | 409 | 23.3 | 2.78 | 4.92 | 3.76 | 3.60 | 3.44 | 3.29
${}_{\;\,42}^{100}\mathrm{Mo}$ | $\rightarrow$ | ${}_{\;\,44}^{100}\mathrm{Ru}$ | 3034.40(17) | [46] | 0.564 | 428 | 241 | 419 | 24.1 | 2.94 | 5.17 | 4.00 | 3.82 | 3.62 | 3.46
${}_{\;\,48}^{116}\mathrm{Cd}$ | $\rightarrow$ | ${}_{\;\,50}^{116}\mathrm{Sn}$ | 2813.50(13) | [47] | 0.601 | 451 | 229 | 442 | 22.0 | 3.42 | 5.92 | 4.74 | 4.51 | 4.17 | 3.97
${}_{\;\,52}^{128}\mathrm{Te}\,$ | $\rightarrow$ | ${}_{\;\,54}^{128}\mathrm{Xe}$ | 865.87(131) | [48] | 0.589 | 452 | 206 | 457 | 33.1 | 3.74 | 6.42 | 5.29 | 5.04 | 4.53 | 4.32
${}_{\;\,52}^{130}\mathrm{Te}\,$ | $\rightarrow$ | ${}_{\;\,54}^{130}\mathrm{Xe}$ | 2526.97(23) | [47] | 0.589 | 452 | 206 | 457 | 33.1 | 3.74 | 6.42 | 5.29 | 5.04 | 4.53 | 4.32
${}_{\;\,54}^{136}\mathrm{Xe}$ | $\rightarrow$ | ${}_{\;\,56}^{136}\mathrm{Ba}$ | 2457.83(37) | [49] | 0.606 | 476 | 217 | 465 | 15.2 | 3.91 | 6.67 | 5.57 | 5.31 | 4.71 | 4.49
${}_{\;\,60}^{150}\mathrm{Nd}$ | $\rightarrow$ | ${}_{\;\,62}^{150}\mathrm{Sm}$ | 3371.38(20) | [50] | 0.519 |  |  | 514 | 16.7 |  |  | 6.50 | 6.20 |  |
${}_{\;\,92}^{238}\mathrm{U}\;\;$ | $\rightarrow$ | ${}_{\;\,94}^{238}\mathrm{Pu}$ | 1437.3 | [51] | 0.546 |  |  | 774 | 17.5 |  |  | 14.58 | 13.90 |  |

The Hamiltonian of the daughter ion’s electrons is related to the Hamiltonian
of the parent neutral atom’s electrons via the relation

$\hat{H}_{Z+2,Z}=\hat{H}_{Z,Z}-2\sum_{i}\frac{1}{{r}_{i}},$ (1)

where $r_{i}=|\mathbf{r}_{i}|$, $\mathbf{r}_{i}$ is the coordinate of the
$i$th electron, and summation is performed by $i=1,\ldots,Z$. The electrons of
the daughter ion are in the state $|Z,Z\rangle$ for the next moment after the
decay, while the nucleus acquires a charge of $Z+2$. The relationship

$\mathcal{M}=\langle Z,Z|\hat{H}_{Z+2,Z}|Z,Z\rangle-\langle
Z+2,Z|\hat{H}_{Z+2,Z}|Z+2,Z\rangle$ (2)

determines the average excitation energy of the daughter ion’s electrons, or,
with account of Eq. (1),

$\mathcal{M}=E_{Z,Z}+2Z^{-1}E_{Z,Z}^{\mathrm{C}}-E_{Z+2,Z},$ (3)

where $E_{Z,N}^{\mathrm{C}}$ is the Coulomb interaction energy of the
electrons with the nucleus.

Table 1 shows the results of the calculation of the excitation energy in the
TF, TFDW and DHF models. First, the values $\mathcal{M}^{\prime}$ are found,
which differ from $\mathcal{M}$ by replacing in Eq. (3) the binding energy of
the electrons of the ion $E_{Z+2,Z}$ with the binding energy of the electrons
of the neutral atom $E_{Z+2,Z+2}$. The difference between $E_{Z+2,Z}$ and
$E_{Z+2,Z+2}$ is equal to the double ionization energy, $I_{2}$; there is a
relation $\mathcal{M}=\mathcal{M}^{\prime}-I_{2}$. The experimental values of
$I_{2}$ are collected in [38].

In the TF model, the calculations are carried out according to the scheme of
[15]. The TFWD model, being a generalization of the TF model, additionally
takes into account exchange contribution to the energy of electron gas [25]
and spatial inhomogeneity in the electron density [26]. A consistent
semiclassical decomposition of the density functional with account of the
inhomogeneity can be found in the monograph by Kirzhnits [27]. In its simplest
form

$\displaystyle E_{Z,N}$ $\displaystyle=$ $\displaystyle\int
d\mathbf{r}\left(-\frac{Z}{r}n(\mathbf{r})+c_{1}n^{5/3}(\mathbf{r})+c_{2}n^{4/3}(\mathbf{r})+c_{3}\frac{(\nabla
n(\mathbf{r}))^{2}}{n(\mathbf{r})}\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{1}{2}\int
d\mathbf{r}d\mathbf{r}^{\prime}n(\mathbf{r})\frac{1}{|\mathbf{r}-\mathbf{r}^{\prime}|}n(\mathbf{r}^{\prime}).$

Here, $n(\mathbf{r})$ is the electron density, the first term under the
integral sign represents the interaction energy of the electrons with the
nuclus, $E_{Z,N}^{\textrm{C}}$, the second term is the kinetic energy, the
third one is the exchange energy, the fourth one is the Weizsacker gradient
correction [26]. The last term is the interaction energy of electrons. The
coefficients $c_{i}$ equal

$c_{1}=\frac{3}{10}(3\pi^{2})^{2/3},~{}~{}c_{2}=-\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3},~{}~{}c_{3}=\frac{\lambda}{8}.$
(5)

The value $\lambda=1/5$ of the phenomenological models [28, 29] is used.

In [28], the binding energy of neutral atoms N, Ne, Ar, Kr, Xe with filled
valence shells is calculated using the TFDW model. Parameterization of the
results gives $E_{Z,Z}=-0.536Z^{2.38}$, which is not much different from the
TF model, where $E_{Z,Z}=-0.764Z^{7/3}$. The energy of the Coulomb interaction
of electrons with the nucleus is calculated using the screening function.
Integrating the expression for $E_{Z,Z}^{\textrm{C}}$ by parts, the action of
the Laplacian is transferred to the Coulomb potential, which gives a delta
function at the origin. The difference between the total potential and the
nuclear potential occurs as a multiplier. The interaction energy turns out to
be $Z^{2}(a-b)$, the screening function parameters $a$ and $b$ are given in
Table II of [28]. The fitting gives $E_{Z,Z}^{\textrm{C}}=-1.270Z^{2.38}$ in
agreement with the virial theorem. The parameterization accuracy is not worse
than 0.5%. The corresponding results for $\mathcal{M}$ are shown in Table 1.

The average values of $\langle Z,Z|r^{-1}_{i}|Z,Z\rangle$ required to estimate
$E_{Z,Z}^{\mathrm{C}}$ in the DHF method are tabulated in [30, 31, 32, 33]. In
[33], the values of $E_{Z,Z}^{\mathrm{C}}$ are also provided. The results of
calculations of $\mathcal{M}$ within the framework of DHF model [33] are shown
in Table 1.

The TF and DHF models agree well with each other and agree qualitatively with
the predictions of the TFDW model.

The variance of the electron excitation energy is determined by the formula

$\mathcal{D}=\langle Z,Z|\hat{H}_{Z+2}^{2}|Z,Z\rangle-\langle
Z,Z|\hat{H}_{Z+2}|Z,Z\rangle^{2}.$ (6)

Taking into account Eq. (1) we have

$\frac{1}{4}\mathcal{D}=\sum_{ij}\langle
Z,Z|\frac{1}{r_{i}}\frac{1}{r_{j}}|Z,Z\rangle-\langle
Z,Z|\sum_{i}\frac{1}{r_{i}}|Z,Z\rangle^{2}.$ (7)

The summation is performed in the range $1\leq i,j\leq Z$. In the TF and TFDW
models, the two-particle electron density is not defined, however, it is
possible to fix the upper limit of the variance [15]:

$\frac{1}{4}\bar{\mathcal{D}}=\int
d\mathbf{r}\frac{1}{{r}^{2}}n(\mathbf{r})-Z^{-1}\left(\int
d\mathbf{r}\frac{1}{{r}}n(\mathbf{r})\right)^{2}.$ (8)

Calculation of the integral of $1/r^{2}$ over the electron density
distribution in the TFDW model leads to values that can be parameterized as

$\int d\mathbf{r}\frac{1}{r^{2}}n(\mathbf{r})=5.81Z^{2.00}.$ (9)

The parameterization accuracy is not worse than 5%. The values of
$\bar{\mathcal{D}}$ in the TF and TFDW models are shown in Table 1.

In the DHF method, it is possible to estimate not only the upper bound of the
variance, but also the variance itself. In disregard of exchange effects
$\mathcal{D}$ is calculated from Eq. (7) after factorization of the average
value under the double summation sign. The corresponding results, using the
tabulated values of averages $1/r_{i}$ and $1/r^{2}_{i}$ for the electron
orbitals [32], are shown in Table 1.

The exchange effects are taken into account by averaging the two-particle
operator over the total wave function of electrons of the atom. In the one-
determinant approximation, the wave function has the form

$\Psi_{\alpha_{1}\alpha_{2}\ldots\alpha_{N}}=\frac{1}{\sqrt{N!}}\epsilon^{s_{1}s_{2}...s_{N}}\phi_{\alpha_{s_{1}}}^{1}\phi_{\alpha_{s_{2}}}^{2}\ldots\phi_{\alpha_{s_{N}}}^{N},$
(10)

where $\phi_{\alpha}^{i}$ are the wave functions of electrons, the index
$i=1,...,N$ counts the spatial coordinates and spin indices, the index
$\alpha$ counts the quantum numbers of orbitals. In the case under
consideration, $\alpha=(njlm)$, where $n$ is the principal quantum number, $j$
is the total angular momentum, $m$ is its projection, $l=j\pm 1/2$ is the
orbital angular momentum. A fixed set of quantum numbers
$(\alpha_{1},\alpha_{2},\ldots,\alpha_{N})$ determines the state of the
electron shells of the atom. The tensor $\epsilon^{s_{1}s_{2}\ldots s_{N}}=\pm
1$ performs antisymmetrization.

The functions $\phi_{\alpha}^{i}$ are orthonormal. We write them as the
product of the radial and angular parts:

$\phi_{njlm}^{i}=R_{njl}(r_{i})\Omega_{jm}^{l}(\mathbf{n}_{i}).$ (11)

Here $R_{njl}(r)$ is a real function, $\Omega_{jm}^{l}(\mathbf{n})$ is a
spherical spinor depending on the unit vector
$\mathbf{n}=\mathbf{r}/|\mathbf{r}|$. We denote by $\kappa_{njl}$ the number
of occupied energy levels with quantum numbers $(njl)$. In the case of fully
occupied energy levels, as well as cases allowing for each pair of $(jl)$ the
existence of no more than one partially occupied energy level with the maximum
total angular momentum, $j^{\max}=\kappa_{njl}(2j+1-\kappa_{njl})/2$, one can
simplify Eq. (7) by replacing the summation over electrons by the summation
over energy levels:

$\frac{1}{4}\mathcal{D}=\sum_{njl}\kappa_{njl}\langle
njl|r^{-2}|njl\rangle-\sum_{nn^{\prime}jl}\min(\kappa_{njl},\kappa_{n^{\prime}jl})\langle
njl|r^{-1}|n^{\prime}jl\rangle^{2}.$ (12)

The matrix elements are defined according to

$\langle njl|h(r)|n^{\prime}jl\rangle=\int
r^{2}drh(r)R_{{}_{njl}}(r)R_{{}_{n^{\prime}jl}}(r).$ (13)

The sum of the diagonal components $n=n^{\prime}$ of Eq. (12) coincides with
the right side of Eq. (7) through factorization of the mean value under the
sign of the double sum, as it is assumed in the TF and TFDW estimates. The
components with $n\neq n^{\prime}$ in the second term of Eq. (12) are related
to exchange effects. The exchange effects reduce the variance.

In the RHF method, the functions $R_{njl}(r)$ are tabulated [30]. To calculate
the variance taking into account exchange effects, knowledge of the off-
diagonal matrix elements $\langle njl|r^{-1}|n^{\prime}jl\rangle$ and $\langle
njl|r^{-2}|n^{\prime}jl\rangle$ is required. Table 2 shows the results for
molybdenum atom in the RHF method. The diagonal matrix elements are compared
with those in the DHF method [32]. There is some systematic underestimation of
the average values in comparison with the DHF method, which is due to the
shift in relativistic models of the electron density to smaller distances
[35]. A similar pattern holds for other 10 atoms. Accordingly, the variance in
the RHF method without taking into account exchange contributions is also
systematically lower than the predictions of the DHF method.

The average values of $1/r$ and $1/r^{2}$ in inner and outer orbitals are
approximately in the ratios $Z:1$ and $Z^{2}:1$, which is consistent with the
values of the diagonal matrix elements in Table 2. In cases where for a
partially occupied level the total angular momentum is not the maximum and/or
where there exist more than one partially occupied level for a pair of $(jl)$,
the formula (12) is used as an approximation. Since in medium-heavy and heavy
atoms, the main contribution to the variance comes from electrons in inner
shells, where $\kappa_{njl}=2j+1$, one can expect that accuracy of such an
estimate is quite high.

The results of calculations of $\mathcal{D}$ in the RHF method, taking into
account exchange contributions, are presented in Table 1. For comparison, the
results of calculations without exchange effect are also provided. The
agreement with the TF, TFDW, RHF and DHF models is quite satisfactory.

For applications, we recommend the values of excitation energy
$\mathcal{M}_{\mathrm{DHF}}$, for variance – $\mathcal{D}_{\mathrm{DHF/b}}$,
as theoretically the most justified. The estimate of
$\mathcal{D}_{\mathrm{DHF/b}}$ differs from $\mathcal{D}_{\mathrm{DHF/a}}$ in
that it includes exchange corrections found by the RHF method. Taking into
account various approximations, the uncertainty in
$\mathcal{M}_{\mathrm{DHF}}$ and $\mathcal{D}_{\mathrm{DHF/b}}$ can be
estimated at $<10\%$.

In the non-relativistic TF, TFDW and RHF models $\mathcal{M}$ only weakly
depends on $Z$, while $\mathcal{D}$ grows approximately as $Z^{2}$. This
behavior is perfectly consistent with the highlighted role of K electrons,
whose nonrelativistic excitation theory in $\beta$ decay is developed in [39,
40].

Table 2: Matrix elements $\langle nl|r^{-1}|n^{\prime}jl\rangle$ and $\langle
njl|r^{-2}|n^{\prime}jl\rangle$ for $n\leq n^{\prime}$ electron orbitals in a
molybdenum atom. Calculations use electron wave functions of the RHF method
[30] with degeneracy in $j$. In the lower part of table, diagonal matrix
elements $n=n^{\prime}$ of the relativistic DHF method [32] are given; the
upper and lower rows of P and D waves correspond to $j=l+1/2$ and $j=l-1/2$,
respectively.  42Mo
---
| $\langle njl|r^{-1}|n^{\prime}jl\rangle$ | 1S | 2S | 3S | 4S | 5S |  | 2P | 3P | 4P |  | 3D | 4D
---|---|---|---|---|---|---|---|---|---|---|---|---
1S | 41.49 | 7.962 | 3.231 | -1.255 | 0.321 |  |  |  |  |  |  |
2S |  | 9.378 | 2.160 | -0.803 | 0.204 |  2P | 9.339 | -1.858 | -0.626 |  |  |
3S |  |  | 3.264 | -0.665 | 0.163 |  3P |  | 3.164 | 0.582 |  3D | 2.970 | -0.361
4S |  |  |  | 1.171 | -0.149 |  4P |  |  | 1.052 |  4D |  | 0.714
5S |  |  |  |  | 0.327 |  |  |  |  |  |  |
[32] | 43.55 | 9.939 | 3.409 | 1.209 | 0.322 |  | 9.412 | 3.190 | 1.059 |  | 2.958 | 0.695
|  |  |  |  |  |  | 9.879 | 3.300 | 1.089 |  | 2.987 | 0.705
$\langle njl|r^{-2}|n^{\prime}jl\rangle$ | 1S | 2S | 3S | 4S | 5S |  | 2P | 3P | 4P |  | 3D | 4D
1S | 3455. | 984.9 | 410.3 | -160.1 | 40.94 |  |  |  |  |  |  |
2S |  | 357.8 | 141.7 | -55.02 | 14.06 |  2P | 118.4 | -37.69 | -13.17 |  |  |
3S |  |  | 65.20 | -24.42 | 6.223 |  3P |  | 21.34 | 6.697 |  3D | 11.17 | -2.120
4S |  |  |  | 10.41 | -2.564 |  4P |  |  | 3.157 |  4D |  | 0.965
5S |  |  |  |  | 0.748 |  |  |  |  |  |  |
[32] | 4005. | 439.4 | 80.03 | 12.73 | 0.830 |  | 120.7 | 21.93 | 3.243 |  | 11.11 | 0.930
|  |  |  |  |  |  | 141.5 | 25.50 | 3.744 |  | 11.37 | 0.960

The parameter $K_{Z}$ given in Table 1 represents the overlap amplitude of the
wave functions of all the electrons in the ground state of the parent atom
with the wave functions of the electrons in the ground state of the twice
ionized daughter atom, whose electrons have retained their initial
configuration. The corresponding wave functions of the electrons are not
orthogonal because the charges of the nuclei before and after shaking differ
by two units, and as a result, the overlap of electron wave functions with the
identical quantum numbers is not equal to one. The daughter ion gets excited
as a result. The value $K^{2}_{Z}$ determines the probability of inheriting
quantum numbers by the electrons and, accordingly, the absence of shaking
effects.

To estimate $K_{Z}$, a multiparticle calculation using the DHF method is
required, which was performed using the software package Grasp-2018 [34, 35].
A set of large $f_{njl}^{+}(r)$ and small $f_{njl}^{-}(r)$ radial components
of electron wave functions for all quantum numbers $(njl)$ is obtained for
each parent atom of Table 1 with an appropriate electron configuration and a
total angular momentum corresponding to the ground state of the parent atom’s
electrons. Similarly, for a daughter ion with a nuclear charge $Z+2$, a set of
radial components $\tilde{f}_{njl}^{\pm}(r)$ is obtained. The overlap
amplitude $\mathcal{O}_{njl}$ of the wave functions of electrons with the same
quantum numbers is equal to

$\mathcal{O}_{njl}=\int\left(\tilde{f}_{njl}^{+}(r)f_{njl}^{+}(r)+\tilde{f}_{njl}^{-}(r)f_{njl}^{-}(r)\right)r^{2}dr.$

In the one-determinant approximation and without taking into account exchange
terms, the amplitude of $K_{Z}$ is equal to the product of the amplitudes of
$\mathcal{O}_{njl}$ in the degree equal to the occupation number of the
corresponding level:

$K_{Z}=\prod_{njl}\left(\mathcal{O}_{njl}\right)^{\kappa_{njl}}.$ (14)

For a wide range of atomic numbers $Z$, the values of $K_{Z}$ turn out to be
close to 1/2. The probability of $K^{2}_{Z}\sim 1/4$ is quite small, which
indicates the dominance of channels with the excited electron shells of atoms
in agreement with the phenomenological analysis [20].

The above approach assumes that the resulting configuration of the daughter
ion is the ground state with quantum numbers of electrons inherited from the
parent atom. However spectroscopic analysis shows that this condition is not
always met. For example, the Ti III ion formed after $0\nu 2\beta$ decay of Ca
with the electron configuration [Ar]4s2, in the ground state has the
configuration [Ar]3d2. Strictly speaking, this fact means that the overlap is
exactly zero: $K_{Z}\equiv 0$, therefore, the decay with the dominant
probability is accompanied by the excitation of the electron shells of the
atom. The energy of the [Ar]4s2 lowest configuration exceeds that of the
[Ar]3d2 configuration of Ti III by 12.7 eV. A similar situation occurs in
double-$\beta$ decay atoms of Zr, Mo, Nd and U. In these cases, the amplitude
$K_{Z}$, given in Table 1, is the amplitude of the transition to the most
likely excited state of the electron shells of the daughter ion. In
approximately every fourth case, double-$\beta$ decay of Ca, Zr, Mo, Nd and U
is accompanied by the de-excitation of the electron shells of atoms from the
unique excited state to the ground state with the emission of a series of
photons of the ultraviolet range. The observation of these photons, whose
wavelengths are well known, can serve as an auxiliary signature for the
identification of decay.

Knowledge of the parameters $K_{Z}$, $\mathcal{M}$ and $\mathcal{D}$ is
sufficient to construct simple models of the energy distribution of
$\beta$-electrons in $0\nu 2\beta$ decay. With a probability of $K_{Z}^{2}$,
the electrons of the decaying atom remain in the lowest energy state,
preserving their quantum numbers, with a probability of $1-K_{Z}^{2}$ they
pass into an excited state. The conditional probability of transition to an
excited state with energy $\epsilon$ in the interval $d\epsilon$ is denoted by
$w(\epsilon/Q^{*})d\epsilon/Q^{*}$. The total probability density takes the
form

$p(\epsilon)=K_{Z}^{2}\delta(\epsilon)+(1-K_{Z}^{2})w(\epsilon/Q^{*})/Q^{*}.$
(15)

The binomial distribution is used for $w(x)$, which has a certain
universatility and is widely used in modeling random processes [41]. The
distribution has two free parameters, which are fixed by normalization to the
average value of $\mathcal{M}=\int_{0}^{Q^{*}}d\epsilon\epsilon p(\epsilon)$
and the mean square of the energy
$\mathcal{D}+\mathcal{M}^{2}=\int_{0}^{Q^{*}}d\epsilon\epsilon^{2}p(\epsilon)$.

Based on the DHF model predictions, we calculate the maximum deviation,
$\Delta T_{\max}$, of the $\beta$-electrons energy from the decay energy
$Q^{*}$. $\Delta T_{\max}$ can be determined through the equation

$p_{T}=\int_{0}^{\Delta T_{\max}}d\epsilon p(\epsilon)$

for a given probability, $p_{T}$. The value of $p_{T}=0.9$ corresponds to the
deviations of the $\beta$-electrons energy from $Q^{*}$ less than $\Delta
T_{\max}=180$ eV (Ca), 18 eV (Ge), 19 eV (Se), and $\Delta T_{\max}<5$ eV for
Zr, Mo, Cd, Te, Xe, Nd, U. At the probability of $p_{T}=0.95$, the deviations
do not exceed $\Delta T_{\max}=1.16$ keV (Ca), 0.55 keV (Ge), 0.44 keV (Se),
0.25 keV (Zr), 0.18 keV (Mo), 69 eV (Cd), 22 eV (128Te), 30 eV (130Te), 19 eV
(Xe), 11 eV (Nd), and $\Delta T_{\max}<5$ eV for U.

The decay energy without the energy taken away by neutrinos is measured in
calorimetric detectors, where the energy resolution reaches a few keV (GERDA).
Using a track calorimeter, the SuperNEMO experiment measures the energy of
$\beta$-electrons in $0\nu 2\beta$-selenium decay with an uncertainty of 4% at
an energy of $Q$. Innovative technologies with excellent energy resolution are
in high demand for reducing background noise and for tracking the impact of
atomic shell excitations on the neutrino mass constraints.

To summarize, the overlap amplitudes for the electron shells of the parent
atom and the daughter ion for each atom whose $2\nu 2\beta$ decay was observed
experimentally were found. In the double-$\beta$ decay of atoms 82Se, 96Zr,
100Mo, 150Nd, and 238U, the electron shells with probability $\sim 1/4$ turn
out to be the lowest excited state with quantum numbers inherited from the
parent atom. Such decays are accompanied by a subsequent de-excitation with
characteristic emission of photons of the ultraviolet range. In the atoms
48Ca, 76Ge, 116Cd, 128Te, 130Te, and 136Xe, the daughter ion’s electrons move
to the ground state with a probability of $\sim 1/4$ and to an excited state
with a probability of $\sim 3/4$. The average value and variance of the
excitation energy were computed for each of the scenarios under consideration.
The dependence on the atomic number indicates the dominant contribution to the
variance of the Feinberg–Migdal effect. Deviations of the $\beta$-electrons
energy from the decay energy $Q^{*}$ were estimated for the neutrinoless mode
of double-$\beta$ decay.

The work was supported by the grant #23-22-00307 of Russian Science
Foundation.

## References

  * [1] D. G. Phillips II et al., Phys. Rep. 612, 1–45 (2016).
  * [2] J. Schechter and J. W. F. Valle, Phys. Rev. D 25, 2951 (1982).
  * [3] M. Hirsch, S. Kovalenko, and I. Schmidt, Phys. Lett. B 642, 106–110 (2006).
  * [4] S. Weinberg, Phys. Rev. Lett. 43, 1566 (1979).
  * [5] The GERDA Collaboration, Nature 544, 47 - 52 (2017).
  * [6] G. Anton et al. [EXO-200 Collaboration], Phys. Rev. Lett. 123, 161802 (2019).
  * [7] R. Arnold et al. [NEMO-3 Collaboration], Phys. Rev. D 89, 111101(R) (2014).
  * [8] R. Arnold et al., [NEMO-3 Collaboration], Phys. Rev. D 92, 072011 (2015)
  * [9] D. Q. Adams et al. [CUORE Collaboration], Phys. Rev. Lett. 124, 122501 (2020).
  * [10] D. Q. Adams, C. Alduino, K. Alfonso, F. T. Avignone, O. Azzolini, G. Bari, F. Bellini, G. Benato, M. Beretta, M. Biassoni, A. Branca [CUORE Collaboration], Nature 604, 53–58 (2022).
  * [11] S. Abe, S. Asami, M. Eizuka, S. Futagi, A. Gando, Y. Gando, T. Gima, A. Goto, T. Hachiya, K. Hata, S. Hayashida, K. Hosokawa, K. Ichimura, S. Ieki [KamLAND-Zen Collaboration], Phys. Rev. Lett. 130, 051801 (2023).
  * [12] F. Šimkovic, Amand Faessler, V. Rodin, P. Vogel, J. Engel, Phys. Rev. C 77, 045503 (2008).
  * [13] J. T. Suhonen, Front. Phys. 5, 55 (2017).
  * [14] F. Šimkovic, Phys. Usp. 64, 1238 (2021).
  * [15] M. I. Krivoruchenkoa, K. S. Tyrin, and F. F. Karpeshin, Pis’ma v ZhETF 117, 887–891 (2023) [JETP Lett. 117, 884–888 (2023)].
  * [16] I. Lindgren, Journal of Electron Spectroscopy and Related Phenomena 137-140, 59–71 (2004).
  * [17] D. S. Akerib et al. (LUX Collaboration), Phys. Rev. Lett. 122, 131301 (2019).
  * [18] E. Aprile et al. (XENON Collaboration), Phys. Rev. Lett. 123, 241803 (2019).
  * [19] P. Agnes et al., (DarkSide Collaboration), Phys. Rev. Lett. 130, 101001 (2023).
  * [20] M. I. Krivoruchenko, K. S. Tyrin, Eur. Phys. J. A 56, 16 (2020).
  * [21] F. F. Karpeshin, M. B. Trzhaskovskaya, and L. F. Vitushkin, Yad. Fiz. 83, 344 (2020). [Phys. At. Nucl. 83, 608 (2020)].
  * [22] F. F. Karpeshin, M. B. Trzhaskovskaya, Yad. Fiz. 85, 387 (2022). [Phys. At. Nucl. 85 (2020)].
  * [23] F. F. Karpeshin and M. B. Trzhaskovskaya, Phys. Rev. C 107, 045502 (2023).
  * [24] L. D. Landau, E. M. Lifschitz, Quantum Mechanics: Non-relativistic Theory. Course of Theoretical Physics. Vol. 3 (3rd ed.), (London: Pergamon1977).
  * [25] P. A. M. Dirac, Math. Proc. Cambridge Phil. Soc. 26, 376–385 (1930).
  * [26] C. F. v. Weizsäcker, Zeitschrift für Physik 96, 431 (1935).
  * [27] D. A. Kirzhnits, Field Theoretical Methods in Many-Body Systems, (Pergamon Press, Oxford 1967), – pp. 394.
  * [28] E. K. U. Gross and R. M. Dreizler, Phys. Rev. A 20, 1798 (1979).
  * [29] W. Stich, E. K. U. Gross, P. Malzacher, and R.M. Dreizler, Z. Phys. A 309, 5-11 (1982).
  * [30] E. Clementi and C. Roetti, At. Data Nucl. Data Tables 14, 177–478 (1974).
  * [31] C. C. Lu, T. A. Carlson, F. B. Malik, T. C. Tucker, and C. W. Nestor, Jr., At. Data Nucl. Data Tables 3, 1 (1971).
  * [32] J. P. Desclaux, At. Data Nucl. Data Tables, 12, 31l–406 (1973).
  * [33] Ken-Ning Huang, Michio Aoyagi, Mau Hsiung Chen, Bern Grasemann, Hans Mark, At. Data Nucl. Data Tables 18, 243–291 (1976).
  * [34] K. G. Dyall, I. P. Grant, C. T. Johnson, F. A. Parpia, and E. P. Blummer, Comput. Phys. Commun. 55, 425–456 (1989).
  * [35] I. P. Grant, Relativistic Quantum Theory of Atoms and Molecules: Theory and Computation, (Springer Science+Business Media, New York 2007).
  * [36] I. M. Band, M. A. Listengarten, M. B. Trzhaskovskaya, and V. I. Fomichev, Computer Program Complex RAINE I–IV, Leningrad Nuclear Physics Institute Reports LNPI-289 (1976), LNPI-298 (1977), LNPI-299 (1977), and LNPI-300 (1977).
  * [37] I. M. Band, M. B. Trzhaskovskaya, C. W. Nestor Jr., P. O. Tikkanen, and S. Raman, At. Data Nucl. Data Tables 81, 1 (2002).
  * [38] A. Kramida, Yu. Ralchenko, J. Reader and NIST ASD Team (2022). NIST Atomic Spectra Database (ver. 5.10), https://physics.nist.gov/asd. National Institute of Standards and Technology, Gaithersburg, MD. DOI: https://doi.org/10.18434/T4W30F
  * [39] E. L. Feinberg, J. Phys. (USSR) 4, 423 (1941).
  * [40] A. Migdal, ZhETF 11, 207–212 (1941).
  * [41] V. S. Korolyuk, N. I. Portenko, A. V. Skorokhod, A. F. Turbin, Handbook of Probability Theory and Mathematical Statistics (Moscow, Nauka 1985). – 640 pp. [in Russian].
  * [42] A. A. Kwiatkowski, T. Brunner, J. D. Holt, A. Chaudhuri, U. Chowdhury, M. Eibach, J. Engel, A. T. Gallant, A. Grossheim, M. Horoi, A. Lennarz, T. D. Macdonald, M. R. Pearson, B. E. Schultz, M. C. Simon, R.A. Senkov, V. V. Simon, K. Zuber, J. Dilling, Phys. Rev. C 89, 045502 (2014).
  * [43] M. Suhonen, I. Bergström, T. Fritioff, Sz. Nagy, A. Solders and R. Schuch, J. Instrum. 2 06003, (2007).
  * [44] David L. Lincoln, Jason D. Holt, Georg Bollen, Maxime Brodeur, Scott Bustabad, Jonathan Engel, Samuel J. Novario, Matthew Redshaw, Ryan Ringle, and Stefan Schwarz, Phys. Rev. Lett. 110, 012501 (2013).
  * [45] M. Alanssari, D. Frekers, T. Eronen, L. Canete, J. Dilling, M. Haaranen, J. Hakala, M. Holl, M. Jeskovsky, A. Jokinen, A. Kankainen, J. Koponen, A. J. Mayer, I. Moore, D. Nesterenko, I. Pohjalainen, P. Povinec, J. Reinikainen, S. Rinta-Antila, P. C. Srivastava, Jouni Suhonen, R. I. Thompson, Annika Voss, M. E. Wieser, Phys. Rev. Lett. 116, 072501 (2016).
  * [46] S. Rahaman, V.-V. Elomaa, T. Eronen, J. Hakala, A. Jokinen, J. Julin, A. Kankainen, A. Saastamoinen, J. Suhonen, C. Weber, J. Aysto, Phys. Lett. B 662, 111–116 (2008)
  * [47] S. Rahaman, V.-V. Elomaa, T. Eronen, J. Hakala, A. Jokinen, A. Kankainen, J. Rissanen, J. Suhonen, C. Weber, J. Aysto, Phys. Lett. B 703, 412–416 (2011).
  * [48] N. D. Scielzo, S. Caldwell, G. Savard, J. A. Clark, C. M. Deibel, J. Fallis, S. Gulick, D. Lascar, A. F. Levand, G. Li, J. Mintz, E. B. Norman, K. S. Sharma, M. Sternberg, T. Sun, J. Van Schelt, Phys. Rev. C 80, 025501 (2009).
  * [49] Matthew Redshaw, Elizabeth Wingfield, Joseph McDaniel, and Edmund G. Myers, Phys. Rev. Lett. 98, 053003 (2007).
  * [50] V. S. Kolhinen, T. Eronen, D. Gorelov, J. Hakala, A. Jokinen, A. Kankainen, I. D. Moore, J. Rissanen, A. Saastamoinen, J. Suhonen, and J. Äystö, Phys. Rev. C 82, 022501 (2010).
  * [51] Richard B. Firestone, Virginia S. Shirley, editors; Coral M. Baglin, S.Y. Frank Chu, and Jean Zipkin, assistant editors. Table of isotopes (Wiley-Interscience, New York 1996).


above process on the new sub-path until we find $X(v)$. To make the strategy
possible it suffices to build a binary search data structure for every
majority path.

For any node $v^{\prime}$ of the aggregation-tree $T$, let $T(v^{\prime})$
denote the subtree of $T$ rooted at $v^{\prime}$. For every majority path
$P_{m}$, we build a binary tree $T_{P_{m}}$ for its nodes, say
${v_{1},v_{2},\ldots,v_{m}}$. First we assign a weight to each node $v_{i}$.
Let $\\{v_{1}^{\prime},v_{2}^{\prime}\ldots,v_{k}^{\prime}\\}$ be all children
of $v_{i}$ that are not in the majority path $P_{m}$. The weight of $v_{i}$ is
$1$ plus the total size of
$T(v_{1}^{\prime}),T(v_{2}^{\prime}),\ldots,T(v_{k}^{\prime})$, where the size
of a subtree is the number of its nodes. To build the binary tree $T_{P_{m}}$
for ${v_{1},v_{2},\ldots,v_{m}}$, we first find the weighted median node
$v_{j^{\prime}}$ and make $v_{j^{\prime}}$ the root of $T_{P_{m}}$; then
recursively build a subtree for ${v_{1},v_{2},\ldots,v_{j^{\prime}-1}}$ and
let its root be the left child of $v_{j^{\prime}}$; also recursively build a
subtree for ${v_{j^{\prime}+1},v_{j^{\prime}+2},\ldots,v_{m}}$ and let its
root be the right child of $v_{j^{\prime}}$.

Let $Z^{\prime}(v^{\prime})$ be the sub-path of $Z(v)$ starting at node
$v^{\prime}\in Z(v)$ and ending at the last node $X(v)$ of $Z(v)$. Consider
the following FindTail procedure.

Algorithm 7 FindTail($T_{s}$, $v_{s}$)

Input: A subtree $T_{s}$ of the aggregation-tree $T$ with its root being the
head of a majority path in $T$. A node $v_{s}$ of $T$ in $Z(v)$.
Output: $X(v)$.

1:Let $T_{P_{m}(v_{s})}$ be the binary tree built for the majority path
$P_{m}(v_{s})$ containing $v_{s}$. Conduct a binary search on
$T_{P_{m}(v_{s})}$ to find the last node $w$ in the majority path
$P_{m}(v_{s})$ which also appears in $Z^{\prime}(v_{s})$.

2:If all children of $w$ in $T$ are not in $Z(v)$, return $w$ as the last node
$X(v)$ of $Z(v)$.

3:Otherwise, exactly one child of $w$, say $w_{Z}$, is in $Z(v)$. Call
FindTail($T(w_{Z})$, $w_{Z}$).

(Note: Input $T_{s}$ above is for purpose of analysis. When FindTail($T_{s}$,
$v_{s}$) is called, it means $X(v)$ is found to be in $T_{s}$ and FindTail
will search in $T_{s}$ for $X(v)$. However $T_{s}$ is not actually used in the
procedure.)

Call FindTail($T$, $v$) to find the last node $X(v)$ of $Z(v)$.

As stated earlier, for any descendant node $v^{\prime}$ of $v$ in the
aggregation-tree $T$, it takes $O(1)$ time to decide whether $v^{\prime}\in
Z(v)$. Using this as a basic decision operation, in the procedure FindTail
above, the binary search in Step 1 is performed as follows. Start at the root
$v_{r,s}$ of the binary tree $T_{P_{m}(v_{s})}$. If $v_{r,s}$ is a $Z(v)$ node
and is also the last $Z(v)$ node in the majority path $P_{m}(v_{s})$, then we
are done with the binary search on $T_{P_{m}(v_{s})}$. If $v_{r,s}$ is a
$Z(v)$ node but is not the last $Z(v)$ node in the majority path
$P_{m}(v_{s})$ (i.e., this can be decided by checking whether the child
$v^{\prime}$ of $v_{r,s}$ in $T$ along the path $P_{m}(v_{s})$ appears in
$Z(v)$), then search recursively on the right child of $v_{r,s}$ in the binary
tree $T_{P_{m}(v_{s})}$. If $v_{r,s}$ is not a $Z(v)$ node, then search
recursively on the left child of $v_{r,s}$ in the binary tree. From this, it
is clear that calling FindTail($T$, $v$) will eventually find the last node
$X(v)$ of $Z(v)$.

Clearly, in the procedure FindTail($T_{s}$, $v_{s}$), $T(w_{Z})$ is at most of
half the size of $T_{s}$ due to the property of a majority path. Thus, after
each recursion of FindTail, the search space is reduced by at least half of
the size. Since the aggregation-tree $T$ has $O(n)$ nodes, FindTail takes
$O(\log n)$ recursions.

Let $W$ be the weight of the node $w$ ($w$ is found in Step 1 of
FindTail($T_{s}$, $v_{s}$)) and $W^{\prime}$ be the size of $T_{s}$. Clearly,
$W^{\prime}$ is equal to the total weight of nodes of $T_{P_{m}(v_{s})}$
(since the root of $T_{s}$ is the head of $P_{m}(v_{s})$). It is easy to see
that the binary search in Step 1 takes $O(1)+O(\log\frac{W^{\prime}}{W})$
time. Also, $W$ is larger than the size of $T(w_{Z})$. FindTail($T$, $v$)
produces a sequence of recursive calls. Let $W_{1},W_{2},\ldots,W_{a}$ and
$W_{1}^{\prime},W_{2}^{\prime},\ldots,W_{a}^{\prime}$ be the values of $W$ and
$W^{\prime}$, respectively, in the sequence of FindTail calls, sorted by the
time of the calls. Note that $a=O(\log n)$ by the above discussion. Then the
running time of FindTail($T$, $v$) is

$a\times
O(1)+O(\log\frac{W^{\prime}_{1}}{W_{1}}+\cdots+\log\frac{W^{\prime}_{a}}{W_{a}})$
$=O(\log n)+O(\log W^{\prime}_{1}-\log W_{1}+\log W^{\prime}_{2}-\log
W_{2}+\cdots+\log W^{\prime}_{a}-\log W_{a})\leq O(\log n)+O(\log
W^{\prime}_{1}).$

The last inequality follows from the fact that $W^{\prime}_{i+1}\leq W_{i}$.
Since $W_{1}$ is the size of $T$, which is $O(n)$, it then follows that the
running time of FindTail($T$, $v$) is $O(\log n)$.

The time and space of the preprocessed data structure are clearly $O(n\log n)$
and $O(n)$, respectively. Thus the lemma is true. ∎

###### Lemma 13

If $Y(v)$ is not empty, then $v_{l}$ is either $X(v)$ or a descendant of
$X(v)$ in $T$.

###### Proof

First, we show that $v_{l}$ is either a node in $Z(v)$ or a descendant of
$X(v)$. Suppose this is not the case. Let $v_{P}$ be the last node of $Z(v)$
such that $v_{P}$ is an ancestor of $v_{l}$ in $T$. Since $v\in Z(v)$ is an
ancestor of $v_{l}$, such a node $v_{P}$ must exist. Since $v_{P}$ is not
$X(v)$ and $v_{P}$ is not $v_{l}$, there are two distinct children of $v_{P}$,
say $v_{1}$ and $v_{2}$, such that $v_{1}$ is $v_{l}$ or an ancestor of
$v_{l}$, and $v_{2}$ is $X(v)$ or an ancestor of $X(v)$. This means that both
$R^{\prime}(v_{1})$ and $R^{\prime}(v_{2})$ intersect $B$ (by the definitions
of $v_{l}$ and $X(v)$). But this contradicts with the fact that $v_{P}$ is
non-splittable.

Next, we show that if $R(v_{l})$ is disjoint with $B$ or $S(v_{l})\leq\lVert
q_{c}-L(v_{l})\rVert\Delta^{-1}(\epsilon)/(3d)$, then $R(v_{l})$ does not
fully contain all edges in $Y(v)$. The case where $R(v_{l})$ is disjoint with
$B$ is trivial. Thus we focus only on the case of $S(v_{l})\leq\lVert
q_{c}-L(v_{l})\rVert\Delta^{-1}(\epsilon)/(3d)$. Suppose by contradiction
$R(v_{l})$ fully contains all edge in $Y(v)$. Then the edge length of
$R(v_{l})$ is larger than that of $B$. Recall that $\Delta^{-1}(\epsilon)$ is
set to be no bigger than $1/2$. It is impossible that the edge length of
$R(v_{l})$ is larger than that of $B$ and also satisfies the inequality
$S(v_{l})\leq\lVert q_{c}-L(v_{l})\rVert\Delta^{-1}(\epsilon)/(3d)$. The
reason is the following. Since $R(v_{l})$ intersects $B$, $\lVert
q_{c}-L(v_{l})\rVert$ is no larger than the edge length $D(B)$ of $B$ plus the
diameter of $R(v_{l})$. From this, we know that $S(v_{l})\leq\lVert
q_{c}-L(v_{l})\rVert\Delta^{-1}(\epsilon)/(3d)$ implies

$S(v_{l})\leq(D(B)+dS(v_{l}))\Delta^{-1}(\epsilon)/(3d)\leq(S(v_{l})+dS(v_{l}))\Delta^{-1}(\epsilon)/(3d)\leq(1+d)S(v_{l})/(6d).$

This is impossible for any $d\geq 1$. Therefore, if $R(v_{l})$ is disjoint
with $B$ or $S(v_{l})\leq\lVert
q_{c}-l(v_{l})\rVert\Delta^{-1}(\epsilon)/(3d)$, then $v_{l}$ is not in $Z(v)$
and must be a descendant of $X(v)$.

Finally, if $R(v_{l})$ intersects $B$ and $S(v_{l})>\lVert
q_{c}-l(v_{l})\rVert\Delta^{-1}(\epsilon)/(3d)$, then $v_{l}$ must be
splittable (otherwise, $v_{l}$ will not be the last node in Step 5 of
SlowFind). Suppose $v_{l}$ is not $X(v)$ or a descendant of $X(v)$. Then
$v_{l}$ and one of its children must be in $Z(v)$. This means that $S(v_{l})$
is at least 2 times the edge length of $B$. Since $v_{l}$ is splittable, when
$R(v_{l})$ is divided into $2^{d}$ sub-boxes (in a quad-tree decomposition),
at least one of these sub-boxes has one facet, say $f$, which intersects $B$
and is inside $R(v_{l})$ (i.e. is not part of a face of $R(v_{l})$) . The
facet $f$ must intersect one of $B$’s edges, because its edge length is no
smaller than that of $B$. Therefore, some edge $e$ of $B$ must be cut after
the decomposition. Note that $e$ cannot be any edge in $Y(v)$, since otherwise
no child of $v_{l}$ will fully contain $e$, and this contradicts with the fact
that one child of $v_{l}$ is in $Z(v)$ and fully contains all edges in $Y(v)$.
This also means that $e$ is not entirely in $R(v_{l})$. Therefore, $e$ must
pass through one of the $2^{d}$ sub-boxes of $R(v_{l})$ so that it can be
possibly cut. This implies that the length of $e$ is larger than half of
$S(v_{l})$, which is a contradiction. Hence, $v_{l}$ is $X(v)$ or a descendant
of $X(v)$. ∎

The above lemmas suggest that if $Y(v)$ is not empty, then we can use FindTail
to first find $X(v)$ of $Z(v)$. $v_{l}$ is either $X(v)$ or its descendant. If
it is the first case, then we have already found $v_{l}$. Otherwise, it means
that at least one of the edges in $Y(v)$ has been cut while decomposing the
box of $X(v)$. Thus, we can first determine the child $v^{\prime}$ of $X(v)$
which is $v_{l}$ or its ancestor (i.e., using the fact that
$R^{\prime}(v^{\prime})$ intersects $B$). Then we generate a new set of edges,
$Y(v^{\prime})$, of $B$ which are fully contained in $R(v^{\prime})$. Clearly,
the size of $Y(v^{\prime})$ is reduced by at least $1$ from that of $Y(v)$. If
$Y(v^{\prime})$ is not empty, then we repeat the above procedure to find a new
$X(v^{\prime})$. Since the size of $Y(v)$ is a constant, after a constant
number of iterations, it will become zero. At that time, we can use the binary
search method on $P_{I}(v^{\prime\prime})$ to eventually find $v_{l}$, where
$v^{\prime\prime}$ is the last node from the above process. The total time of
the entire process is only $O(\log n)$. This leads us to the following
improved procedure Find($v,q_{c}$) for finding an effective cover.

Algorithm 8 Find($v$, $q_{c}$)

Input: A node $v$ in the aggregation-tree $T$ and a query point $q_{c}$.
Output: Part of an effective cover for $q_{c}$ in the subtree $T(v)$ of $T$.

1:If $R(v)$ does not intersect $B$, return.

2:If $R(v)$ is small enough, i.e., $S(v)\leq\lVert
q_{c}-L(v)\rVert\Delta^{-1}(\epsilon)/(3d)$, report $v$ as one of the output
nodes, return.

3:If $v$ is splittable, call Find($v_{i}$, $q_{c}$) on each of $v$’s children,
$v_{i}$, in $T$ that touches $B$, return.

4:Let $v_{0}=v$ and $i=0$. While $Y(v_{i})$ is not empty, do

  1. a.

Use FindTail to find $X(v_{i})$. If $X(v_{i})$ is $v_{l}$, let $v=X(v_{i})$
and go to Step 1.

  2. b.

Otherwise, let $v_{i+1}$ be the child of $X(v_{i})$ which is $v_{l}$ or its
ancestor. Let $i=i+1$ and continue the while loop.

Use binary search on $P_{I}(v_{i})$ to find $v_{l}$. Let $v=v_{l}$ and go to
Step 1.

###  5.5 Algorithm Analysis

Now we analyze the running time of our assignment algorithm.

###### Lemma 14

Step 4 of Find$(v,q_{c})$ takes $O(\log n)$ time.

###### Proof

Based on the above discussions, we know that the while loop in Step 4 can
execute at most $O(1)$ iterations. In each iteration, the time is dominated by
that of FindTail, which takes $O(\log n)$ time. Thus, the total time of the
while loop is $O(\log n)$. The binary search on $P_{I}(v_{i})$ takes $O(\log
n)$ time. Hence, the lemma follows. ∎

The next lemma bounds the total time of the procedure Find($v,q_{c}$).

###### Lemma 15

An effective cover of size $O(\log n)$ for $q_{c}$ can be obtained by the
procedure Find($v_{r},q_{c}$) in $O(\log^{2}n)$ time.

###### Proof

Since Algorithm Find improves only Step 5 of SlowFind, its recursion tree is
the same as SlowFind and thus is of size $O(\log n)$ (by Lemma 11). By Lemma
14, we know that each recursion of Find takes $O(\log n)$ time. Hence, the
total time for finding an effective cover is $O(\log^{2}n)$. ∎

###### Lemma 16

Algorithm Assign$(c)$ takes $O(\log^{\max\\{2,d\\}}n)$ time to assign a
$(1-\epsilon)$-approximate maximum influence site to each type-2 cell of the
AI decomposition.

###### Proof

The running time of Step 2 is $O(\log^{2}n)$ (by Lemma 15). For the running
time of Step 3, since an effective cover is of size $O(\log n)$ and each
partition induced by a hyperplane can be determined by $q_{c}$ and $d-1$ nodes
(or more precisely, their representative points) in the cover, the total
number of these partitions is hence $O(\log^{d-1}n)$. Each partition takes
$O(\log n)$ time to compute the influence. Thus the total time of Step 3 is
$O(\log^{\max\\{2,d\\}}n)$. Other steps take $O(1)$ time. Therefore, the total
time of Algorithm Assign is $O(\log^{\max\\{2,d\\}}n)$. The quality guarantee
follows from Lemma 10. ∎

###### Theorem 5.2

A $(1-\epsilon)$-approximate vector CIVD can be constructed in
$O(n\log^{\max\\{3,d+1\\}}n)$ time, where $n$ is the number of input points in
$P$ and $d$ is the dimensionality of the space.

###### Proof

The correctness and approximation ratio follow from Lemma 16. For the running
time, we know that the AI decomposition takes $O(n\log n)$ time to generate
totally $O(n\log n)$ cells. Each cell takes $O(\log^{\max\\{2,d\\}}n)$ time to
determine its approximate maximum influence site. Other preprocessing takes
$O(n\log n)$ time. Thus, the total time is $O(n\log^{\max\\{3,d+1\\}}n)$. ∎

##  6 Density-based CIVD

In this section, we show how to augment the AI-Decomposition algorithm to
generate a $(1-\epsilon)$-approximate CIVD for the density-based CIVD problem.

###  6.1 Problem Description and Properties of the Influence Function

The density-based CIVD problem for a set $P$ of $n$ points in $\mathbb{R}^{d}$
is to partition the space into cells so that all points in each cell share the
same subset $C$ of $P$ as their densest cluster (see Fig. 17). For a given
query point $q\in\mathbb{R}^{d}$, the densest cluster $C_{m}(P,q)$ of $q$ is
the subset $C$ of $P$ which maximizes the influence $F(C,q)=\lvert
C\rvert/V(C,q)$ over all subsets of $P$, where
$V(C,q)=\frac{\pi^{\frac{d}{2}}l^{d}}{\Gamma(\frac{d}{2}+1)}$ is the volume of
the smallest ball centered at $q$ and containing all points in $C$, $l$ is the
maximum distance from $q$ to any point in $C$, and $\Gamma$ is the gamma
function in the volume computation of a $d$-dimensional ball. In other words,
$C_{m}(P,q)$ is the cluster with the highest density around $q$. Fig. 18 shows
an example of an approximate density-based CIVD generated by AI Decomposition.

Clearly, density-based CIVD is closely related to the widely used density-
based clustering problem [12, 13, 15, 30, 35]. In the density-based clustering
problem, two parameters (the radius $r$ of the neighborhood ball $B$ and the
density $d$ of the input points inside $B$) are used to partition the input
points into non-overlapping clusters. Different from the density-based
clustering problem, our problem does not use such parameters, and can
automatically determine the radius of each dense cluster. It generates not
only the dense clusters but also their associated Voronoi cells. Since
density-based clustering is used in many data mining, pattern recognition,
biomedical imaging, and social network applications, we expect that the
density-based CIVD is also applicable in these areas. Also, since our problem
allows the generated clusters to overlap with one another, it has the
potential to be applicable to overlapping clustering problems [1, 7, 10, 17].

Figure 17: Examples of the densest clusters for three query points
$q_{1},q_{2}$, and $q_{3}$. Figure 18: An example of an approximate density-
based CIVD for $4$ input points on the plane (generated by our algorithm with
$\beta=0.1$). The figure on the right is a zoomed view of the figure on the
left.

Similar to the vector CIVD problem, our goal for the density-based CIVD is
also a $(1-\epsilon)$-approximate CIVD. To use the AI decomposition for this
problem, we first show that it satisfies the three properties in Section 3.

###### Theorem 6.1

The density-based CIVD problem satisfies the three properties in Section 3.

###### Proof

First, we show Property 2. Consider a set $C$ of input points and a query
point $q$ in $\mathbb{R}^{d}$. Let $\psi$ be an $\epsilon$-perturbation on $P$
(with the witness point $q$) for some constant $0<\epsilon<1$. Let $p_{max}$
be the point in $C$ that is farthest from $q$, and $\psi(p_{max}^{\prime})$ be
the point in $\psi(C)$ that is farthest from $q$. Since

$\lVert\psi(p_{max}^{\prime})-q\rVert\leq(1+\epsilon)\lVert
p_{max}^{\prime}-q\rVert$

and

$\lVert p_{max}^{\prime}-q\rVert\leq\lVert p_{max}-q\rVert,$

we have

$\lVert\psi(p_{max}^{\prime})-q\rVert\leq(1+\epsilon)\lVert p_{max}-q\rVert.$
(6)

Furthermore, from

$\lVert\psi(p_{max}^{\prime})-q\rVert\geq\lVert\psi(p_{max})-q\rVert$

and

$\lVert\psi(p_{max})-q\rVert\geq(1-\epsilon)\lVert p_{max}-q\rVert,$

we get

$\lVert\psi(p_{max}^{\prime})-q\rVert\geq(1-\epsilon)\lVert p_{max}-q\rVert.$
(7)

By the influence function, and inequalities (6) and (7), we know that

$(1+\epsilon)^{-d}F(C,q)\leq F(\psi(C),q)\leq(1-\epsilon)^{-d}F(C,q).$

This implies Property 2 (by setting
$\delta(\epsilon)=\max\\{1-(1+\epsilon)^{-d},(1-\epsilon)^{-d}-1\\}$).

For Property 3, we assume that, for any query point $q\in\mathbb{R}^{d}$,
there is a point $p\in P$ and a subset $A$ of points in $P$ such that for
every $a\in A$, $\lVert a-q\rVert>n^{1/d}\lVert q-p\rVert$. For any subset
$B\subseteq P$ intersecting with $A$, let $b$ be a point in $A\cap B$. Then
$\lVert b-q\rVert>n^{1/d}\lVert q-p\rVert$. Now we compare $F(B,q)$ with
$F(\\{p\\},q)$. It is clear that the smallest ball centered at $q$ and
containing $B$ is at least $n$ times larger (in volume) than the smallest ball
centered at $q$ and containing $\\{p\\}$. Since $\lvert B\rvert\leq n$, we
have $F(B,q)<F(\\{p\\},q)$.

If there is a subset $P^{\prime}\subseteq P$ and $p\in P^{\prime}$ such that
$n^{1/d}\lVert q-p\rVert<\epsilon^{\prime}\cdot\lVert
q-p^{\prime}\rVert<\lVert q-p^{\prime}\rVert$ for all $p^{\prime}\in
P\setminus P^{\prime}$ and some constant $0<\epsilon^{\prime}<1$, then by the
above discussion, we have $C_{m}(P,q)=C_{m}(P^{\prime},q)$. This means that
for every cluster $C$ and query point $q$, the pair $(C,q)$ is stable. Thus
Property 3 holds.

For Property 1, it is clear that after a scaling or a rotation about any query
point $q\in\mathbb{R}^{d}$, the distance from every point in $P$ to $q$ is
changed by the same factor which is uniquely determined by the transformation.
From the influence function, we know that Property 1 holds. ∎

###  6.2 Assignment Algorithm by Modifying the AI Decomposition

To make use of the AI decomposition to construct an approximate density-based
CIVD, our idea is to modify the AI-Decomposition algorithm (Algorithm 3) so
that some additional information is maintained for assigning a cluster to each
resulted type-2 cell. (Note that by Theorem 4.1, for each type-1 cell $c$, we
can simply use the distance-node $v$ which dominates $c$ as its densest
cluster.) In this way, we can obtain the approximate CIVD at the same time
when completing the AI decomposition.

Recall that an input point $p$ is recorded in the AI-Decomposition algorithm
only when its distance to the current to-be-decomposed box is large enough.
Therefore, for a cell $c$, it is most likely that an input point recorded
earlier is farther away from $c$ than an input point recorded later.
Intuitively, the recorded distances (of the input points) should be roughly in
a decreasing order with respect to the order in which they are recorded. Below
we discuss how to utilize this observation to modify the AI-Decomposition
algorithm. A proof of this observation will be given later.

To show how to modify the AI-Decomposition algorithm, we first consider an
example. Let $q$ be a query point, and $P=\\{p_{1},p_{2},\ldots,p_{n}\\}$ be a
set of input points in the decreasing order of their distances to $q$ (i.e.,
$\lVert p_{i}-q\rVert>\lVert p_{j}-q\rVert$ for all $1\leq i<j\leq n$, and no
two different points in $P$ have the same distance to $q$). To find
$C_{m}(P,q)$, we can use the following approach which scans $P$ only once in
its sorted order and uses $O(1)$ additional space. For each $1\leq i\leq n$,
we compute $D_{i}=\frac{c_{d}(n-i+1)}{\lVert p_{i}-q\rVert^{d}}$, and store
the largest $D_{i}$ during the scanning process, where
$c_{d}=\frac{\Gamma(\frac{d}{2}+1)}{\pi^{\frac{d}{2}}}$. Since the ball
centered at $q$ and with radius $\lVert p_{i}-q\rVert$ contains exactly
$n-i+1$ input points, $\\{p_{i},p_{i+1},\ldots,p_{n}\\}$, the largest $D_{i}$
value, along with the corresponding $i$, gives us the desired densest cluster
$C_{m}(P,q)=\\{p_{i},p_{i+1},\ldots,p_{n}\\}$.

Figure 19: A configuration with 22 input points, $B(u)$ being processed,
$v_{1}$ and $v_{2}$ been removed from $L$ before, and $v$ being removed from
$L$. For any $q$ in $B$ if we draw a ball centered at $q$ with radius $r$
being roughly the distance between $q$ and $v$($\lVert q-l(v)\rVert$), the
ball should “approximately” include all the input points not yet been
recorded.(i.e. points not in $v_{1}$ and $v_{2}$.) Denote set of these points
by $P_{u}$. Even without knowing much information about $P_{u}$ at this point,
it is possible to obtain an approximate value of $F(P_{u},q)$.
($\frac{c_{2}(22-3-4)}{r^{2}}$ in this case.) Since this works for all $q$ in
$B(u)$, This information can be passed down during the recursion.

With the above illustration, we can now modify the Decomposition algorithm
(Algorithm 2), as follows. In particular, we change Step 2 of the
Decomposition algorithm, since this is the step in which distance-nodes are
removed from the list $L$. Before the execution of Step 2, we sort the
distance-nodes in $L$ by the decreasing order of their distances to the
current box-node $u$ (if multiple nodes have the same distance, then we order
them arbitrarily). Then, we execute Step 2 and try to remove distance-nodes
from $L$ according to this sorted order. We assume that in the Decomposition
algorithm, a number $M$ is maintained for storing the total number of input
points which are recorded for the current box-node $u$. During the execution
of Step 2, after removing each distance-node $v$, we compute a value
$D=\frac{c_{d}(n-M)}{r^{d}}$ and then update $M$ (i.e., increase $M$ by the
cardinality $|P_{v}|$ of $v$). We save the largest value $D$ along with the
corresponding distance-node $v$ and the box-node $u$. In each recursive call
to the Decomposition algorithm, we pass the stored $D$, $u$, and $v$ to the
next level of the recursion.

See Figure 19 for better understanding of the strategy.

Clearly, such a modification on Step 2 of the Decomposition algorithm
resembles the computation in the above example. The only difference is that in
the above example, the input points are considered strictly in the decreasing
order of their distances to the query point $q$, but in the modified
Decomposition algorithm, distance-nodes are not always removed by the
decreasing order of their distances to some query point $q$. This is because
at different recursion levels, distance-nodes may not be removed in a strictly
decreasing order. Below we show that an approximate densest cluster for $q$
can still be obtained, despite the above difference.

Let $c$ be a type-2 cell generated by the AI decomposition and $q$ be a query
point in $c$. Consider the root-to-$c$ path in the recursion tree of the
Decomposition algorithm for $c$. Let $v_{1},v_{2},\ldots,v_{m}$ be the
sequence of distance-nodes removed from $L$ along this recursion path (sorted
by the increasing order of the time when they are removed), and
$x_{1},x_{2},\ldots,x_{m}$ be the closest distances to their corresponding
box-nodes $u$ at the time when they are removed. Let $D_{max}$ be the maximum
value of $D$ passing through this recursion path and $v_{max}$, and $u_{max}$
be the corresponding box-node and distance-node when $D$ achieves its maximum
value.

Below we prove a claim that if $v_{max}=v_{i}$ for some $i$, then the union of
$v_{i},v_{i+1},\ldots,v_{m}$ is almost the densest cluster for $q$ for a
properly chosen $\beta$. Since $v_{1},v_{2},\ldots,v_{m}$ are all distance-
nodes recorded for the type-2 cell $c$, by Lemma 3, we know that they form a
partition of $P$. For any $p\in v_{j}$, by Lemma 2, we have

$(1-\beta)x_{j}\leq\lVert q-p\rVert\leq(1+\beta)x_{j}.$ (8)

Let $\psi$ be a mapping defined as follows: For any $p\in v_{j}$, $\psi(p)$ is
on the ray that emits from $q$ and passes through $p$, and with $\lVert
p-q\rVert=x_{j}$. Let $C$ denote the union of $v_{i},v_{i+1},\ldots,v_{m}$. By
Lemma 1, we know that to prove the above claim, it is sufficient to show that
$F(\psi(C),q)$ is almost as large as $F(C_{m}(P^{\prime},q),q)$, where
$P^{\prime}=\psi(P)$. Clearly, $P^{\prime}$ can be partitioned into subsets
$\psi(v_{1}),\psi(v_{2}),\ldots,\psi(v_{m})$, with all points in each
$\psi(v_{i})$, for $i=1,\ldots,m$, having the same distance $x_{i}$ to $q$.
Based on the above discussion, we know that if $x_{1},x_{2},\ldots,x_{m}$ are
in decreasing order, then $\psi(C)$ is exactly $C_{m}(P^{\prime},q)$. The
following lemma shows that $x_{1},x_{2},\ldots,x_{m}$ are actually in a
roughly sorted order, which is sufficient for us to obtain an approximate
densest cluster.

###### Lemma 17

In the modified AI-Decomposition algorithm with an error tolerance $\beta$,
$x_{j}\leq(1+\beta)x_{i}$ for any $1\leq i<j\leq m$.

###### Proof

From the above discussion, we know that if $v_{i}$ and $v_{j}$ are removed in
the same recursion of the Decomposition algorithm, then $x_{j}\leq x_{i}$
(since in the same recursion, all distance-nodes are removed in a decreasing
order). Thus we can assume that $v_{i}$ and $v_{j}$ are removed in different
recursions with recursive calls Decomposition$(u_{1},\beta,L_{1},T_{p},r_{1})$
and Decomposition$(u_{2},\beta,L_{2},T_{p},r_{2})$, respectively, where
$u_{1}$ is a proper ancestor of $u_{2}$ in the box-tree $T_{q}$. By
definition, we know that $x_{i}$ is the closest distance between $l(v_{i})$
and $B(u_{1})$, and $x_{j}$ is the closest distance between $l(v_{j})$ and
$B(u_{2})$. Also, by the Decomposition algorithm, we know that $B(u_{2})$ is
contained in $B(u_{1})$.

Consider the execution of Decomposition$(u_{1},\beta,L_{1},T_{p},r_{1})$. Let
$v_{j}^{\prime}$ be the node in $L_{1}$ containing $l(v_{j})$. Clearly,
$v_{j}^{\prime}$ is not removed in Step 2, since otherwise $v_{j}$ will not
appear in Decomposition$(u_{2},\beta,L_{2}$, $T_{p},r_{2})$. This means that
$x_{i}$ is larger than the closest distance $x_{j}^{\prime}$ between
$l(v_{j}^{\prime})$ and $B(u_{1})$, since $v_{i}$ is removed from $L_{1}$ in
Step 2 but $v_{j}^{\prime}$ is not. Let $D(u_{1})$ denote the diameter of
$B(u_{1})$, and $x^{\prime\prime}_{j}$ denote the closest distance between
$B(u_{2})$ and $l(v_{j}^{\prime})$. By the triangle inequality, we have
$x_{j}^{\prime\prime}\leq D(u_{1})+x_{j}^{\prime}$. Since $x_{j}^{\prime}\leq
x_{i}$ and $D(u_{1})\leq x_{i}\beta/2$ (by the fact that $v_{i}$ is removed
from $L$ in Step 2), we have $x_{j}^{\prime\prime}\leq(1+\beta/2)x_{i}$.

Let $E^{\prime}(v_{j}^{\prime})$ be the box co-centered with
$E(v_{j}^{\prime})$ and with the edge length half of that of
$E(v_{j}^{\prime})$ (see Algorithm 1). Consider the following two possible
cases.

  1. 1.

_$B(u_{1})$ does not intersect $E^{\prime}(v_{j}^{\prime})$._ In this case, we
have

$\lVert l(v_{j})-l(v_{j}^{\prime})\rVert\leq s(v_{j}^{\prime})\leq
x_{j}^{\prime}\beta/2\leq x_{i}\beta/2.$

Note that

$x_{j}\leq\lVert l(v_{j})-l(v_{j}^{\prime})\rVert+x_{j}^{\prime\prime}.$

Thus, we have

$x_{j}\leq(1+\beta)x_{i}.$

  2. 2.

_$B(u_{1})$ intersects $E^{\prime}(v_{j}^{\prime})$._ In this case, note that
some part of $B(u_{1})$ must be outside of $E(v_{j}^{\prime})$, since
otherwise $v_{j}^{\prime}$ would have been deleted in Step 1. From this, we
know that $D(u_{1})\geq\sqrt{d}\alpha/2$, where $\alpha$ is the edge length of
$E^{\prime}(v_{j}^{\prime})$. Since $\beta<1/2$ (by the assumption on
$\beta$), we have

$\lVert l(v_{j})-l(v_{j}^{\prime})\rVert\leq
s(v_{j}^{\prime})=(\alpha/2)\cdot\beta/2\leq\alpha/8.$

Note that since $B(u_{1})$ intersects $E^{\prime}(v_{j}^{\prime})$, we have
$x_{j}^{\prime}\leq\sqrt{d}\alpha/2$. Thus, $x_{j}^{\prime}+\lVert
l(v_{j})-l(v_{j}^{\prime})\rVert=\sqrt{d}\alpha/2+\alpha/8\leq\sqrt{d}\alpha\leq
2D(u_{1})$. By the triangle inequality, we have $x_{j}^{\prime}+\lVert
l(v_{j})-l(v_{j}^{\prime})\rVert+D(u_{1})\geq x_{j}$, which means $x_{j}\leq
3D(u_{1})$. Recall that $D(u_{1})\leq x_{i}\beta/2\leq x_{i}/4$; we have
$x_{i}\geq 4D(u_{1})$. Therefore $x_{i}>x_{j}$.

∎

With the above lemma, we immediately have the following lemma (in which the
notation was defined before Lemma 17).

###### Lemma 18

$F(\psi(C),q)\geq(1+\beta)^{-d}F(C_{m}(P^{\prime},q),q)$.

###### Proof

From the discussion before Lemma 17, we know that $P^{\prime}$ can be
partitioned into $\psi(v_{1}),\psi(v_{2}),\ldots,\psi(v_{m})$, and for any
$p\in\psi(v_{r})$, $\lVert p-q\rVert=x_{r}$ for all $1\leq r\leq m$.

It is easy to see that $C_{m}(P^{\prime},q)$ can be written as the union of
all distance-nodes in $\\{\psi(v_{j})\mid x_{j}\leq x_{i_{max}}\\}$ for some
$1\leq i_{max}\leq m$. $\psi(C)$ is the union of
$\psi(v_{i}),\psi(v_{i+1}),\ldots,\psi(v_{m})$. Let

$F^{\prime}=\frac{c_{d}(n-\lvert v_{1}\rvert-\lvert v_{2}\rvert-\cdots-\lvert
v_{i-1}\rvert)}{x_{i}^{d}}.$

Let $x_{h}$ be the maximum in $\\{x_{i},x_{i+1},\ldots,x_{m}\\}$. Then,

$F(\psi(C),q)=\frac{c_{d}(n-\lvert v_{1}\rvert-\lvert
v_{2}\rvert-\cdots-\lvert v_{i-1}\rvert)}{x_{h}^{d}}.$

By Lemma 17, we know that $x_{h}\leq(1+\beta)x_{i}$. Thus,

$F(\psi(C),q)\geq(1+\beta)^{-d}F^{\prime}.$

Rewrite $\\{\psi(v_{j})\mid x_{j}\leq x_{i_{max}}\\}$ as
$\\{v_{j_{1}},v_{j_{2}},\ldots,v_{j_{w}}\\}$ with $j_{1}<j_{2}<\cdots<j_{w}$.
We have

$F(C_{m}(P^{\prime},q),q)=\frac{c_{d}(\lvert v_{j_{1}}\rvert+\lvert
v_{j_{2}}\rvert+\cdots+\lvert v_{j_{w}}\rvert)}{x_{i_{max}}^{d}}.$

Let

$G^{\prime}=\frac{c_{d}(\lvert v_{j_{1}}\rvert+\lvert
v_{j_{1}+1}\rvert+\cdots+\lvert
v_{m}\rvert)}{x_{j_{1}}^{d}}=\frac{c_{d}(n-\lvert v_{j_{1}-1}\rvert-\lvert
v_{j_{1}-2}\rvert-\cdots-\lvert v_{1}\rvert)}{x_{j_{1}}^{d}}.$

Note that $\lvert v_{j_{1}}\rvert+\lvert v_{j_{2}}\rvert+\cdots+\lvert
v_{j_{w}}\rvert\leq\lvert v_{j_{1}}\rvert+\lvert
v_{j_{1}+1}\rvert+\cdots+\lvert v_{m}\rvert$. Since $x_{i_{max}}\geq
x_{j_{1}}$, we have

$G^{\prime}\geq F(C_{m}(P^{\prime},q),q).$

Since $G^{\prime}\leq F^{\prime}$, it follows that

$F(\psi(C),q)\geq(1+\beta)^{-d}F^{\prime}\geq(1+\beta)^{-d}G^{\prime}\geq(1+\beta)^{-d}F(C_{m}(P^{\prime},q),q).$

∎

Based on the above two lemmas, we immediately have the following theorem.

###### Theorem 6.2

For any $\beta$ satisfying the conditions
$1-(1+\beta)^{-d}\leq\Delta^{-1}(\epsilon)$ and
$\beta\leq\Delta^{-1}(\epsilon)/3$, the modified AI decomposition algorithm
finds a $(1-\epsilon)$-approximate density-based CIVD in $O(n\log^{2}n)$ time.

###### Proof

By Equation (8) and the definition of $\psi$, we know that for each $p\in
v_{k}$,

$\lVert\psi(p)-p\rVert=\lvert\lVert
p-q\rVert-\lVert\psi(p)-q\rVert\rvert=\lvert\lVert p-q\rVert-
x_{k}\rvert\leq\beta
x_{k}\leq\Delta^{-1}(\epsilon)x_{k}/3\leq\Delta^{-1}(\epsilon)\lVert\psi(p)-q\rVert.$

By Lemma 18, we have

$F(\psi(C),q)\geq(1+\beta)^{-d}F(C_{m}(P^{\prime},q),q)\geq(1-\Delta^{-1}(\epsilon))F(C_{m}(P^{\prime},q),q).$

Consider the perturbation $\psi^{-1}$. If $\beta\leq\Delta^{-1}(\epsilon)/3$,
then $\psi^{-1}$ is also a $(\Delta^{-1}(\epsilon)/3)$-perturbation. From the
proof of Theorem 6.1, we know that $(C,q)$ is stable. By Lemma 9 (note that
Lemma 9 still holds for the density-based CIVD problem), we have
$F(C,q^{\prime})\geq(1-\epsilon)F(C_{m}(P,q^{\prime}),q^{\prime})$ for any
point $q^{\prime}$ in $c$.

For the running time, we note that the additional computation in the modified
AI decomposition algorithm includes sorting the distance-nodes in $L$ and
maintaining the values of $M$, $D$, and the corresponding $u$ and $v$.
Clearly, the additional time is dominated by sorting, which takes
$O(|L|\log|L|)$ time for each recursive call to the Decomposition algorithm.
Since the running time of the original Decomposition algorithm is $O(|L|)$,
and the total running time of the entire AI decomposition is the sum over all
recursions of the Decomposition algorithm, the total time of the modified AI
decomposition thus increases only by an $O(\log n)$ factor. Therefore, the
theorem follows. ∎

## References

  * [1] R. Andersen, D.F. Gleich, and V. Mirrokni, “Overlapping Clusters for Distributed Computation,” Proc. 5th ACM International Conference in Web Search and Data Mining, 2012, pp. 273-282.
  * [2] S. Arya and T. Malamatos, “Linear-Size Approximate Voronoi Diagrams,” Proceedings of the 13th annual ACM-SIAM symposium on Discrete algorithms (SODA’02), pp. 147–155, 2002.
  * [3] S. Arya, T. Malamatos, and D. M. Mount, “Space-Efficient Approximate Voronoi Diagrams,” Proc. 34th ACM Symp. on Theory of Computing (STOC 2002), pp. 721–730, 2002.
  * [4] S. Arya, D. M. Mount, N. S. Netanyahu, R. Silverman, and A. Wu, “An Optimal Algorithm for Approximate Nearest Neighbor Searching,” Journal of the ACM, 45 (1998), pp. 891–923.
  * [5] F. Aurenhammer, “Power Diagrams: Properties, Algorithms and Applications,” SIAM J. on Computing, 16(1)(1987), 78-96.
  * [6] F. Aurenhammer, “Voronoi Diagrams – A Survey of a Fundamental Geometric Data Structure,” ACM Computing Surveys, 23(1991), 345-405.
  * [7] A. Banerjee, C. Krumpelman, S. Basu, R.J. Mooney, and J. Ghosh, “Model-based Overlapping Clustering,” Proc. 11th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2005, pp. 532–537.
  * [8] G. Barequet, M.T. Dickerson, and R.L.S. Drysdale III, “2-Point Site Voronoi Diagrams,” Discrete Applied Mathematics, 122(1-3)(2002), 37-54.
  * [9] G. Barequet, M.T. Dickerson, D. Eppstein, D. Hodorkovsky, and K. Vyatkina, “On 2-Site Voronoi Diagrams under Geometric Distance Functions,” Proc. 8th International Symp. on Voronoi Diagrams in Science and Engineering, 2011, pp. 31-38.
  * [10] F. Bonchi, A. Gionis, and A. Ukkonen, “Overlapping Correlation Clustering,” Proc. IEEE 11th International Conference on Data Mining, 2011, pp. 51-60.
  * [11] P. Callahan and R. Kosaraju, “A Decomposition of Multidimensional Point Sets with Applications to $k$-nearest-neighbors and $n$-body Potential Fields,” JACM, 42(1)(1995), 67-90.
  * [12] F. Cao, M. Ester, W. Qian, and A. Zhou, “Density-based Clustering over an Evolving Data Stream with Noise,” Proceedings of the 6th SIAM International Conference on Data Mining, 2006, pp. 328-339.
  * [13] D.Z. Chen, M.H.M. Smid, and Bin Xu, “Geometric Algorithms for Density-based Data Clustering,” Int. J. Comput. Geometry Appl., 15(3)(2005), 239-260.
  * [14] N. Chen, J. Zhu, F. Sun, and E.P. Xing, “Large-margin Predictive Latent Subspace Learning for Multi-view Data Analysis,” IEEE Transaction on Pattern Analysis and Machine Intelligence, 34(12)(2012), 2365-2378.
  * [15] Y. Chen and L. Tu, “Density-based Clustering for Real-time Stream Data,” Proceedings of the 13th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2007, pp. 133-142.
  * [16] C.M. Christoudias, R. Urtasun, and T. Darrell, “Multi-view Learning in the Presence of View Disagreement,” arXiv:1206.3242, June 2012.
  * [17] G. Cleuziou, L. Martin, and C. Vrain, “Poboc: An Overlapping Clustering Algorithm,” Proc. 16th European Conference on Artificial Intelligence, 2004, pp. 440-444.
  * [18] M.T. Dickerson and D. Eppstein, “Animating a Continuous Family of Two-site Voronoi Diagrams (and a Proof of a Bound on the Number of Regions),” Proc. 25th ACM Symp. Computational Geometry, 2009, pp. 92-93.
  * [19] M.T. Dickerson and M.T. Goodrich, “Two-site Voronoi Diagrams in Geographic Networks”, Proc. 16th ACM SIGSPATIAL International Conf. Advances in Geographic Information Systems, 2008, doi:10.1145/1463434.1463504.
  * [20] H. Ding and J. Xu, “Solving Chromatic Cone Clustering via Minimum Spanning Sphere,” Proc. 38th International Colloquium on Automata, Languages and Programming (ICALP), 2011, pp. 773-784.
  * [21] D. Greene and P. Cunningham, “Multi-view Clustering for Mining Heterogeneous Social Network Data,” Proc. 31st European Conference on Information Retrieval, Workshop on Information Retrieval over Social Networks, LNCS, Vol. 5478, 2009.
  * [22] L. Greengard, The Rapid Evaluation of Potential Fields in Particle Systems. MIT Press, Cambridge (1988).
  * [23] L. Greengard. “The Numerical Solution of the N-body Problem,” Computers in Physics, 4, pp. 142–152 (1990).
  * [24] L. Greengard. “Fast Algorithms for Classical Physics.” Science 265, 909–914 (1994).
  * [25] I. Hanniel and G. Barequet, “On the Triangle-Perimeter Two-site Voronoi Diagram,” Trans. on Computational Science, 9(2010), 54-75.
  * [26] S. Har-Peled, “A Replacement for Voronoi Diagrams of Near Linear Size,” Proc. 42nd Annu. IEEE Sympos. Found. Comput. Sci. (FOCS 2001), pp. 94–103, 2001.
  * [27] Sariel Har-Peled and Nirman Kumar, “Down the Rabbit Hole: Robust Proximity Search and Density Estimation in Sublinear Space.” FOCS 2012: 430-439.
  * [28] D. Hodorkovsky, “2-Point Site Voronoi Diagrams,” M.Sc. Thesis, Technion, Haifa, Israel, 2005.
  * [29] P. Indyk and R. Motwani, “Approximate Nearest Neighbors: Towards Removing the Curse of Dimensionality,” Proc. 30th ACM Symp. on Theory of Computing (STOC 1998), pp. 604–613, 1998.
  * [30] H.-P. Kriegel and M. Pfeifle, “Density-based Clustering of Uncertain Data,” Proceedings of the 11th ACM SIGKDD International Conference on Knowledge Discovery in Data Mining, 2005, pp. 672-677.
  * [31] D.T. Lee and R.L.S. Drysdale, III, “Generalization of Voronoi Diagrams in the Plane,” SIAM J. Comput., 10(1)(1981), 73-87.
  * [32] A.Y. Liu and D.N. Lam, “Using Consensus Clustering for Multi-view Anomaly Detection,” Proc. IEEE CS Security and Privacy Workshops, 2012, pp. 117-125.
  * [33] A. Okabe, B. Boots, K. Sugihara, and S.N. Chiu, Spatial Tessellations: Concepts and Applications of Voronoi Diagrams, 2nd Eds., John Wiley & Sons, 2000.
  * [34] E. Papadopoulou, “The Hausdorff Voronoi Diagram of Point Clusters in the Plane,” Algorithmica, 40(2004), 63-82.
  * [35] J. Sander, M. Ester, H.-P. Kriegel, and X. Xu, “Density-Based Clustering in Spatial Databases: The Algorithm GDBSCAN and Its Applications,” Data Mining and Knowledge Discovery, 2(2)(1998), 169-194.
  * [36] K. Vyatkina and G. Barequet, “On 2-Site Voronoi Diagrams under Arithmetic Combinations of Point-to-Point Distances,” Proc. 7th International Symp. Voronoi Diagrams in Science and Engineering, 2010, pp. 33-41.