# Quantum liquids and droplets with low-energy interactions in one dimension

Ivan Morera  Departament de Física Quàntica i Astrofísica, Facultat de Física,
Universitat de Barcelona, E-08028 Barcelona, Spain  Institut de Ciències del
Cosmos, Universitat de Barcelona, ICCUB, Martí i Franquès 1, E-08028
Barcelona, Spain  Bruno Juliá-Díaz  Departament de Física Quàntica i
Astrofísica, Facultat de Física, Universitat de Barcelona, E-08028 Barcelona,
Spain  Institut de Ciències del Cosmos, Universitat de Barcelona, ICCUB, Martí
i Franquès 1, E-08028 Barcelona, Spain  Manuel Valiente
<EMAIL_ADDRESS>Institute for Advanced Study, Tsinghua University,
Beijing 100084, China

###### Abstract

We consider interacting one-dimensional bosons in the universal low-energy
regime. The interactions consist of a combination of attractive and repulsive
parts that can stabilize quantum gases, droplets and liquids. In particular,
we study the role of effective three-body repulsion, in systems with weak
attractive pairwise interactions. Its low-energy description is often argued
to be equivalent to a model including only two-body interactions with non-zero
range. Here, we show that, at zero temperature, the equations of state in both
theories agree quantitatively at low densities for overall repulsion, in the
gas phase, as can be inferred from the $S$-matrix formulation of statistical
mechanics. However, this agreement is absent in the attractive regime, where
universality only occurs in the long-distance properties of quantum droplets.
We develop analytical tools to investigate the properties of the theory, and
obtain astounding agreement with exact numerical calculations using the
density-matrix renormalization group.

Introduction. Recent advances in the preparation, manipulation and observation
of ultracold atomic droplets and liquids Cabrera2018 ; Schmitt2016 ;
Ferrier2016 ; Ferioli2019 ; Bottcher2019 ; Semeghini2018  have sparked renewed
interest in the physics of these systems Cui2018 ; Hu2020a ; Hu2020b ;
Ferioli2020 ; Astrakharchik2018a ; Parisi2020 ; Wang2020 ; Bottcher2021 ;
Morera2020 ; Morera2021 ; Cikojevic2018 ; Cikojevic2020a ; Cikojevic2020b .
This is more so due to their low densities and temperatures, allowing
universal low-energy descriptions that are independent of the short-distance
details of the relevant interactions Bloch2008 . The possibility to
effectively confine these systems to one spatial dimension Kinoshita2004 ;
Paredes2004 , where interaction effects are enhanced Giamarchi , makes
ultracold atoms a promising platform to realize highly controllable strongly
interacting droplets and liquids. With traditional quantum liquids such as
${}^{4}\mathrm{He}$, interatomic potentials that reproduce essentially all of
their experimentally measurable properties are known accurately Kunitski2020 .
In deep contrast, the underlying interactions in ultracold atomic systems are
highly dependent on the particular atomic species, and applied external
fields, such as those involved in magnetic Feshbach resonances Chin2010  and
transversal confinement Olshanii1998 . Hence, it is impractical, if not
impossible, to attempt as accurate a description as in ${}^{4}\mathrm{He}$ for
each realization of an ultracold atomic liquid. Therefore, a universal low-
energy description of these systems, within the effective field theory (EFT)
paradigm Bedaque2002 ; Hammer2013 , is highly desirable.

In the two-body sector, the simplest EFT including scattering length and both
physical and effective ranges has been recently used to describe a one-
dimensional dimerized liquid in an optical lattice Morera2021 . This system is
typically claimed to be equivalent, at low energies, with a theoretically
simpler EFT that includes the two-body scattering length and a three-body
contact interaction Hammer2013 ; Bulgac2002 , the latter being an emergent
property due to three-body processes with two-body interactions that occur off
the energy shell Valiente2019a ; Pricoupenko2020 . In this Letter, we
investigate one-dimensional quantum liquids and droplets at zero temperature
that are described by such low-energy theories. We find that for overall
repulsion the equations of state for the two-body theory and the EFT including
three-body interaction agree well with each other at low densities, in
accordance with the $S$-matrix formulation of statistical mechanics Dashen1969
. However, for overall attraction their droplets and liquid phases are not
equivalent. This fact is proven by studying universal long-distance
asymptotics in both models. We also develop highly non-perturbative,
analytical approximations whose predictions are in excellent agreement with
exact calculations using the density-matrix renormalization group (DMRG).

Low-energy theory. Our first aim is to exactly simulate a one-dimensional
many-body system of bosons with attractive two-body and repulsive three-body
interactions in free space at zero temperature. Ideally, one could apply
ground-state methods, such as ground state quantum Monte Carlo, to the low-
energy Hamiltonian, given by Pastukhov2019 ; Valiente2019a ; Valiente2019b ;
Valiente2020a ; Valiente2020b ; Sekino2018 ; Hou2020 ; Czejdo2020 ; Hou2019 ;
Drut2018 ; Pricoupenko2019 ; Pricoupenko2020 ; Guijarro2018 ; Nishida2018 ;
LPricoupenko2018 ; LPricoupenko2018 ; LPricoupenko2019

$\displaystyle H_{c}$
$\displaystyle=\sum_{i=1}^{N}\frac{p_{i}^{2}}{2m}+g_{0}\sum_{i<j=1}^{N}\delta(x_{i}-x_{j})$
$\displaystyle+g_{3}\sum_{i<j<k=1}^{N}\delta(x_{i}-x_{j})\delta(x_{j}-x_{k}),$
(1)

where $m$ is the mass of the particles, $g_{0}=-2\hbar^{2}/ma$ is the Lieb-
Liniger coupling constant, with $a$ ($>0$) the scattering length, and $g_{3}$
is the bare three-body coupling constant Pastukhov2019 ; Valiente2019a ;
Sekino2018 . With square cutoff ($\Lambda$) regularization, it is given by
$\hbar^{2}/g_{3}(\Lambda)m=\ln|Q_{*}/\sqrt{3}\Lambda|/\sqrt{3}\pi$, with
$Q_{*}$ the three-body momentum scale Valiente2019a  beyond which the low-
energy theory (1) breaks down. If the three-body interaction is repulsive and
$g_{0}=0$, the $T-$matrix exhibits a Landau pole at energy
$E=-\hbar^{2}Q_{*}^{2}/2m$, which prevents the physical $N$-body ($N>2$)
ground state from being explored using ground state methods Sekino2018 ;
Nishida2018 . Below, we solve this issue by discretizing the problem on a
lattice near, but not in the continuum limit, with a genuinely repulsive
three-body force.

The simplest direct discretization of Hamiltonian (1) is given by a
generalized Bose-Hubbard Hamiltonian,

$\displaystyle H$
$\displaystyle=-J\sum_{j}(b_{j+1}^{\dagger}b_{j}+\mathrm{H.c.})+\frac{U_{2}}{2}\sum_{j}n_{j}(n_{j}-1)$
$\displaystyle+\frac{W}{6}\sum_{j}n_{j}(n_{j}-1)(n_{j}-2)+2J\sum_{j}n_{j}.$
(2)

Above, $J=\hbar^{2}/(2md^{2})$ is the hopping strength, with $d$ the lattice
spacing, $U_{2}=g_{0}/d$ is the on-site two-body interaction strength, $W$
($>0$) is a three-body coupling constant to be determined, $b_{j}$
($b_{j}^{\dagger}$) annihilates (creates) a boson at site $j$
($\in\mathbb{Z}$), and $n_{j}=b_{j}^{\dagger}b_{j}$ is the local number
operator. To relate the three-body lattice coupling constant $W$ to the
continuum momentum scale $Q_{*}$, we solve the three-body problem with
$U_{2}=0$ at vanishing total quasimomentum, and match the continuum and
lattice amplitudes at low energies Valiente2019a , obtaining
$W/J=(\beta+\ln|Q_{*}d|/2\sqrt{3})^{-1}$, where $\beta=-0.1956\ldots$ is a
numerical constant of no physical relevance, i.e., it is regularization-
dependent. Hamiltonian (2) can then be used to simulate a continuum repulsive
three-body force at low densities with a finite lattice spacing, provided that
$W/J>0$, i.e., for $|Q_{*}d|>\exp(2\sqrt{3}|\beta|)\approx 8.4$, corresponding
to the weak to moderate coupling regime for the three-body interaction
Pastukhov2019 ; Valiente2019a ; Valiente2019b ; Valiente2020a ; Valiente2020b
; Sekino2018 ; Hou2020 ; Czejdo2020 ; Hou2019 ; Drut2018 ; Pricoupenko2019 ;
Pricoupenko2020 ; Guijarro2018 ; Nishida2018 ; LPricoupenko2018 ;
LPricoupenko2018 ; LPricoupenko2019 .

To test the ability of the lattice Hamiltonian, Eq. (2), to describe the
three-body repulsive side of its continuum counterpart, Eq. (1), we obtain the
zero-temperature equation of state (EoS) of Hamiltonian (2) with $U_{2}=0$ and
$W/J=1.1$ using DMRG, and compare it with the weak-coupling expansion due to
Pastukhov Pastukhov2019

$\displaystyle\frac{E}{N}$
$\displaystyle=\frac{\hbar^{2}\rho^{2}}{3!m}g(\mu)\Bigg{[}1-\frac{4}{\pi}[g(\mu)]^{1/2}$
$\displaystyle+\frac{g(\mu)}{\sqrt{3}\pi}\left\\{\frac{1}{2}\ln\left|\frac{\sqrt{3}\pi}{g(\mu)}\right|-C\right\\}\Bigg{]},$
(3)

where $\mu=\mu(\rho)=\xi/\rho$ is a renormalization scale,
$g(\mu)=-\sqrt{3}\pi/\ln|Q_{*}/\mu|$ is the renormalized coupling constant
Pastukhov2019 , and $C=-9.863\ldots$. The ambiguity, in perturbation theory,
to choose the scale $\mu$ is identical to that in the 2D Bose gas Beane2010 ;
Beane2018 . To obtain predictive power, we match the weak-coupling EoS (3) at
only one value of the density $\rho$ with the EoS at that density obtained
with DMRG. The results are shown in Fig. 1, where astounding agreement is
observed over all ranges of density.

Figure 1: Zero-temperature equation of state for Hamiltonian (2) with
$U_{2}=0$ and $W/J=1.1$ (filled blue dots), weak-coupling EoS, Eq. (3), at
scale $\mu=\exp(\gamma)Q_{*}^{2}/\sqrt{8}\rho$, corresponding to Pastukhov’s
scale Pastukhov2019  (green dashed line), and at the renormalization scale
$\mu=188.464/\rho d^{2}$ (red solid line). Open black squares correspond to
the EoS obtained from the extended Hubbard model at its two-body resonance
(see text).

Quantum liquid. We now turn our attention to the ground state liquid phase of
Hamiltonian (1), and study it both theoretically, with two different
approximations, and numerically using its lattice discretization (2). We begin
by introducing a non-perturbative decoupling approximation for the three-body
interaction of Hamiltonian (1), valid for large particle numbers
($N\to\infty$), while treating the two-body interaction exactly. In this
approximation, the Hamiltonian $H_{D}$ takes the form Supplemental

$H_{D}=\sum_{i=1}^{N}\frac{p_{i}^{2}}{2m}+G_{N}(\Gamma)\sum_{i<j=1}^{N}\delta(x_{i}-x_{j})+C_{N}(\Gamma),$
(4)

which, up to a constant, has the form of the Lieb-Liniger model with coupling
constant $G_{N}(\Gamma)$. In Eq. (4),
$\Gamma=\langle\delta(x_{1}-x_{2})\rangle$,
$G_{N}(\Gamma)=g_{0}+(N-2)g_{3}\Gamma$, and
$C_{N}=-\Gamma^{2}g_{3}N(N-1)(N-2)/3$. In this approximation, $g_{3}$ ceases
to be a bare coupling constant (i.e. cutoff-dependent), and instead becomes a
renormalized one (i.e. finite) to be fixed to reproduce certain physical
properties, such as the equilibrium binding energy. We note that, since the
original Hamiltonian (1) contains a logarithmic anomaly upon renormalization,
the coupling constant $g_{3}$ in this approximation depends on the state. In
particular, in the ground state, it depends on the bulk density of the
droplet, or the particle number. For $N\to\infty$, and in free space
($L=\infty$), that is, at equilibrium, the ground state $\psi_{N}$ and its
energy $E_{0}(N)$ can be obtained exactly using MacGuire’s solution
MacGuire1966 . Solving the problem self-consistently Supplemental , we obtain

$E_{0}(N)\to-\frac{1}{3}\frac{g_{0}^{2}}{g_{3}^{\mathrm{eq}}}N,\hskip
2.84544ptN\to\infty\,,$ (5)

Figure 2: Ratio $e_{0}/g_{0}\rho_{\mathrm{eq}}$ calculated numerically (blue
dots) for Hamiltonian (2) as function of two-body interaction strength. The
three-body interaction strengths for points labelled with $U_{2}/J=-0.195$,
$-0.329$, $-0.499$ and $-0.725$ are given by $W/J=1.28$, $1.4$, $1.48625$ and
$1.508$, respectively. Dashed line ($2/9=0.222\ldots$) is the prediction from
Eq. (8), and dotted line is the standard mean field prediction.

where $g_{3}^{\mathrm{eq}}$ is the value of $g_{3}$ at equilibrium. The
density profile $\rho_{N}(x)$ can be obtained by computing the density with
respect to a fixed center-of-mass coordinate $X$ as Calogero1975

$\rho_{N}(x)=\int\mathrm{d}x_{2}\ldots\mathrm{d}x_{N}\delta(X)|\psi_{N}(x,x_{2},\ldots,x_{N})|^{2}\,.$
(6)

Its bulk density $\rho_{\mathrm{eq}}$, given by
$\rho_{\mathrm{eq}}=\lim_{N\to\infty}\rho_{N}(0)$, is known analytically for
$N\to\infty$ Calogero1975 ; Gertjerenken2012  and, as opposed to the soliton
of the usual Lieb-Liniger model with no three-body repulsion, it is finite
because $G_{N}=O(N^{-2})$ for large $N$. We have

$\rho_{N}(0)=-\frac{1}{4}\frac{mG_{N}}{\hbar^{2}}N(N-1)\to\frac{3}{2}\frac{|g_{0}|}{g_{3}^{\mathrm{eq}}}=\rho_{\mathrm{eq}}\,.$
(7)

Defining $\mathrm{e}_{\mathrm{eq}}=\lim_{N\to\infty}E_{0}/N$, Eqs. (5) and (7)
can be combined to eliminate the (unknown) effective three-body coupling
$g_{3}^{\mathrm{eq}}$ in favour of measurable physical quantities, obtaining

$e_{\mathrm{eq}}=\frac{2}{9}g_{0}\rho_{\mathrm{eq}}\,.$ (8)

The above equation, which is approximate, yet highly non-perturbative, is one
of the main results of this Letter. It predicts a strongly constrained, linear
relation between the equilibrium energy per particle and density for fixed
two-body interaction strength $g_{0}<0$. To test Eq. (8), in Fig. 2, we plot
the relation $e_{\mathrm{eq}}/g_{0}\rho_{\mathrm{eq}}$ calculated for
Hamiltonian (2) using DMRG for a number of different values of the pair
$(U_{2},W)$. The results show that not only is the linear relation in Eq. (8)
a very good approximation, but also the prediction of the proportionality
constant $2/9=0.222\ldots$ is in excellent agreement with the exact results.
Note that mean-field theory, with chemical potential
$\mu=g_{0}\rho+g_{3}\rho^{2}/2$, also predicts a constant value for
$e_{\mathrm{eq}}/|g_{0}|\rho_{\mathrm{eq}}=1/4$, which is $\sim 15-20\%$ off
the numerical values and the prediction of Eq. (8).

The second approximation we introduce is an improved version of mean-field
theory (iMF), by allowing for $g_{3}$ to depend logarithmically on the density
as $g_{3}(\rho)=g_{3}^{\mathrm{eq}}/[\lambda\ln|\rho/\rho_{\mathrm{eq}}|+1]$,
so that it can account for the anomaly non-perturbatively. The value
$g_{3}^{\mathrm{eq}}$ is given by the coupling constant at equilibrium, and is
fixed by $e_{\mathrm{eq}}$, while $\lambda$ is a dimensionless parameter that
is fixed by the equilibrium density $\rho_{\mathrm{eq}}$, given
$e_{\mathrm{eq}}$. To lowest order, integrating the chemical potential gives
for the energy per particle $e_{\mathrm{iMF}}$ in this approximation,

$e_{\mathrm{iMF}}(\rho)\approx\frac{1}{2}g_{0}\rho+\frac{1}{6}g_{3}(\rho)\rho^{2}\,.$
(9)

Given the freedom of scale in logarithmic running of the coupling constant,
the above approximation is sufficient. The relation between the equilibrium
energy and density can also be obtained eliminating $g_{3}^{\mathrm{eq}}$ and
is given by $e_{\mathrm{iMF}}(\rho_{\mathrm{eq}})=\eta
g_{0}\rho_{\mathrm{eq}}$, with $\eta=(1/2)(1-\lambda)/(2-\lambda)$. The
parameter $\lambda$ is universal in the decoupling approximation (see Eq.
(8)), and given by $\lambda=1/5$ ($\eta=2/9$), leaving effectively the iMF
with only one free parameter, $g_{3}$. In Fig. 3, we plot the EoS at zero
temperature for Hamiltonian (2), together with the improved mean-field
approximation, Eq. (9), showing astounding agreement.

Figure 3: Exact zero-temperature equation of state for Hamiltonian (2) with
$U_{2}/J=-0.33$ and $W/J=1.39$ (blue dots), and improved mean-field
approximation, Eq. (9), with $\rho_{\mathrm{eq}}d=0.576$,
$mg_{3}^{\mathrm{eq}}/\hbar^{2}=0.9899$ and $\eta=0.2119$ (see text). Black
open squares are the exact EoS of the extended Hubbard model with $U/J=10$ and
$V/J=-8/5$ (see text).

Universal and non-universal properties. Three-body interactions are always
present in a many-atom system, and can be genuine or emergent from the off-
shell structure of the two-body interactions Valiente2019a ; Pricoupenko2020 ;
Hammer2013 . In general, these are a combination of genuine and emergent. For
repulsive one-dimensional systems, two-body effective range (on-shell) effects
are typically negligible for low densities Astrakharchik2010 . However, the
effects of a non-zero physical range (off-shell) are important for three-body
processes, especially for large scattering lengths Valiente2019a ;
Pricoupenko2020 , which can be rigorously modelled by simple three-body
forces. This can also be explained qualitatively using field redefinitions
Hammer2013  that trade off-shell low-energy interactions in favour of simpler,
on-shell three-particle forces. Since, in this case, the low-energy two- and
three-body amplitudes are essentially identical, the $S$-matrix formulation of
statistical mechanics Dashen1969  implies thermodynamic equivalence in the gas
phase which, at zero temperature, means this is the case for (overall)
repulsive interactions. For attraction, where liquids may be formed, the
results of Ref. Dashen1969  do not apply and, as we shall see, there is no
such equivalence.

To investigate the claims discussed just above, we choose the simplest lattice
Hamiltonian with large scattering length and non-zero effective range, given
by the extended Hubbard model (EHM), corresponding to Hamiltonian (2), with
$U_{2}$ replaced by $U$, $W=0$, plus an interaction term of the form
$\mathcal{V}=V\sum_{j}n_{j}n_{j+1}$. On the two-body resonance ($1/a=0$),
matching the low-energy two-body amplitude for Hamiltonian (2) requires
$U_{2}=0$, while the EHM requires $U/J=-4(V/J)/(2+V/J)$ Valiente2009 ;
Morera2021 . To fit the low-energy three-body amplitude for the EHM, it is
simplest to consider the finite-size ground state energy of Hamiltonian (2)
with three particles, and match the corresponding energy for the EHM
Luscher1986-1 ; Valiente2019a . For the case of Fig. 1, $W/J=1.1$,
corresponding to $U/J=10$ and $V/J=-8/5$ in the EHM. In Fig. 1, we observe
very good agreement between the two EoS for relatively low densities, while
for larger densities two-body on-shell effects appear to dominate in the EHM.

In the attractive case, which features $N$-body bound states for all $N$, it
is possible to fix the locations of the poles of the $S$-matrix (bound state
energies) rather accurately in both models, but not the residues at the poles
without further parametrization. These are related to the so-called asymptotic
normalization coefficient (ANC) $\gamma_{N}$ for a bound state Taylor2000 ;
Luscher1986-2 ; Koenig2012 . For two particles in one dimension, defining the
normalized relative bound state $\psi_{2}(x)$, $\gamma_{2}$ is defined as
$|\gamma_{2}|=\lim_{x_{12}\to\infty}\left[|\psi_{2}(x_{12})|\exp\left(\sqrt{mE_{2}^{(B)}/\hbar^{2}}|x_{12}|\right)\right]$,
where $E_{2}^{(B)}$ ($>0$) is the binding energy. For large and positive
(attractive) scattering length $a$ in comparison with the effective range $r$
($a/|r|\gg 1$), the two-body binding energies obtained with and without
including the effective range agree well with each other. The ANCs are also in
rather good agreement. For example, for the two models considered here, we
have
$\left|\gamma_{2}^{\mathrm{E}}/\gamma_{2}\right|^{2}\approx(r_{e}/d)\left(2-r_{e}/d\right)+\left(2-(r_{e}/d)^{2}\right)/(a/d)$
where $\gamma_{2}^{\mathrm{E}}$ and $\gamma_{2}$ are, respectively, the two-
body ANC of the EHM and Hamiltonian (2) with identical scattering lengths. For
the case studied in Fig. 3, we have
$|\gamma_{2}^{\mathrm{E}}/\gamma_{2}|^{2}\approx 1.07$, which shows a small
yet non-negligible deviation from unity.

For large particle numbers, we can show that the ratios
$|\gamma_{N}^{\mathrm{E}}/\gamma_{N}|$ become exponential in $N$. Using the
universal asymptotics in one dimension Koenig2012 , we obtain Supplemental
for $N$ bosons, as $|x|\to\infty$,

$\rho_{N}(x)\to\frac{\pi}{2}|\gamma_{N}|^{2}A_{N-1}\frac{N}{N-1}e^{-\frac{N}{N-1}2\kappa_{1,N}|x|}.$
(10)

Above, $\gamma_{N}$ is the ANC in the $N\leftrightarrow(N-1)+1$ breakup
channel, $\hbar^{2}\kappa_{1,N}^{2}/2m=|E_{N}-E_{N-1}|(N-1)/N$ represents the
binding energy with respect to the ground state with one fewer particle, and
$A_{N}$ is a model independent normalization factor. Using Eq. (10) for the
two models here considered, we obtain, as $|x|\to\infty$,
$\rho^{\mathrm{E}}_{N}(x)/\rho_{N}(x)\to\left|\gamma_{N}^{\mathrm{E}}/\gamma_{N}\right|^{2}$,
where the superscript E (no superscript) corresponds to the EHM (Hamiltonian
(2)). This relation has important consequences. We assume that the asymptotic
normalization coefficients for the two models are different for all $N$, since
they already are for two particles. If an $N$-body bound state is a quantum
droplet, then most of the particle content lies within the bulk of the
droplet, with a density we may consider constant. The radius $R$ of the
droplet is then $R\approx N/2\rho_{N}(0)$. One may approximate the exponential
tails for large $N$ as $\rho_{N}(x)\sim\rho_{N}(0)\exp(-2\kappa_{1,N}|x-R|)$
as $|x|\to\infty$, obtaining

$\left|\frac{\gamma_{N}^{\mathrm{E}}}{\gamma_{N}}\right|^{2}\sim\frac{\rho_{N}^{\mathrm{E}}(0)}{\rho_{N}(0)}\exp\left[\kappa_{1,N}N\left(\frac{1}{\rho_{N}^{\mathrm{E}}(0)}-\frac{1}{\rho_{N}(0)}\right)\right].$
(11)

The above relation indicates that we require algebraically small
($O(N^{-1/2-|\epsilon|})$) differences between the ANCs in order to achieve
algebraically ($O(N^{-1-|\delta|})$) small differences in the equilibrium
densities of quantum droplets and liquids, and vice versa. A finite difference
in the large-$N$ densities of the two models immediately implies a exponential
disagreement between the ANCs, and therefore the residues at the poles of the
$S$-matrix. Equation (11) shows our claim that, in order for two different
models of a quantum liquid to share thermodynamical properties, not only the
bound state energies but also the residues at these energies must be matched,
and that the matching must have algebraic precision. In Fig. 3, we plot the
EoS for the EHM with $U/J=10$ and $V/J=-8/5$, for which we have verified that
the binding energies for all $N$ match those of Hamiltonian (2) with
$U_{2}=-0.33$ and $W/J=1.39$ within $2\%$. It is observed that the equilibrium
densities are in clear disagreement, and so are the EoS at all densities, as
we expected.

Figure 4: Main panel: Density profiles with $N=50$ particles for Hamiltonian
(2) with $U_{2}/J=-0.33$ and $W/J=1.39$ (open squares) and the extended BH
(see main text) with $U/J=10$ and $V/J=-1.6$ (stars). Inset panel: Density
tails, rescaled by their respective saturation densities. Dashed line
indicates the analytical result Eq. (10).

The disagreement in the EoS between the two models has strong consequences
when looking at the density profiles of the respective quantum droplets. Even
though both models exhibit the same equilibrium energy, their equilibrium
densities, and thus the droplet saturation densities, are different. This
implies that droplets in the two models with equal particle number have
different sizes, see Fig. 4, a direct consequence of the ANCs $\gamma_{N}$,
Eq. (11). On the other hand, the decay of the density far away from the center
of the droplet has to be identical in both models if the equilibrium energy is
the same as predicted by Eq. (10). In the inset panel of Fig. 4 we show that
both droplets exhibit the same exponential density decay far away from the
respective centers. Moreover the decay is dictated by Eq. (10) which shows
that it is directly related to the chemical potential at equilibrium for large
number of particles. Therefore, the tails of quantum droplets are universal
for different models with identical binding energies (in free space), as
opposed to the saturation density.

Conclusions. We have investigated one-dimensional quantum liquids at zero
temperature and droplets using universal low-energy theories. We have shown
that, while in the repulsive case different models with (almost) identical
two- and three-body scattering amplitudes have identical low-density equations
of state, small deviations in densities in the few-particle sector grow
exponentially in the many-body limit. This implies the lack of equivalence of
the zero-temperature equations of state. We have developed theoretical
techniques that yield quantitatively accurate predictions in all density
regimes, when compared to exact results obtained with DMRG.

Acknowledgements This work has been partially supported by MINECO (Spain)
Grant No.FIS2017-87534-P. We acknowledge financial support from Secretaria
d’Universitats i Recerca del Departament d’Empresa i Coneixement de la
Generalitat de Catalunya, co-funded by the European Union Regional Development
Fund within the ERDF Operational Program of Catalunya (project QuantumCat,
ref. 001-P-001644).

## References

  * (1) C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Cheiney and L. Tarruell, Science 359, 301 (2018).
  * (2) M. Schmitt, M. Wenzel, F. Böttcher, I. Ferrier-Barbut and T. Pfau, Nature 539, 259 (2016).
  * (3) I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel and T. Pfau, Phys. Rev. Lett. 116, 215301 (2016).
  * (4) G. Ferioli, G. Semeghini, L. Masi, G. Giusti, G. Modugno, M. Inguscio, A. Gallemí, A. Ricati and M. Fattori, Phys. Rev. Lett. 122, 090401 (2019).
  * (5) G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F. Minardi, M. Modugno, G. Modugno, M. Inguscio and M. Fattori, Phys. Rev. Lett. 120, 235301 (2018).
  * (6) F. Böttcher, J. -N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen and T. Pfau, Phys. Rev. X 9, 011051 (2019).
  * (7) X. Cui, Phys. Rev. A 98, 023630 (2018).
  * (8) H. Hu and X. -J. Liu, Phys. Rev. Lett. 125, 195302 (2020).
  * (9) H. Hu and X. -J. Liu, Phys. Rev. A 102, 043302 (2020).
  * (10) G. Ferioli, G. Semeghini, S. Terradas-Briansó, L. Masi, M. Fattori and M. Modugno, Phys. Rev. Research 2, 013269 (2020).
  * (11) G. E. Astrakharchik and B. A. Malomed, Phys. Rev. A 98, 013631 (2018).
  * (12) L. Parisi and S. Giorgini, Phys. Rev. A 102, 023318 (2020).
  * (13) J. Wang, H. Hu and X. -J. Liu, New J. Phys. 22, 103044 (2020).
  * (14) F. Böttcher, J. -N. Schmidt, J. Hertkorn, K. S. H. Ng, S. D. Graham, M. Guo, T. Langen and T. Pfau, Rep. Progr. Phys. 84, 012403 (2021).
  * (15) I. Morera, G. E. Astrakharchik, A. Polls and B. Juliá-Díaz, Phys. Rev. Research 2, 022008(R) (2020).
  * (16) I. Morera, G. E. Astrakharchik, A. Polls and B. Juliá-Díaz, Phys. Rev. Lett. 126, 023001 (2021).
  * (17) V. Cikojevic, K. Dzelalija, P. Stipanovic, L. V. Markic and J. Boronat, Phys. Rev. B 97, 140502 (2018).
  * (18) V. Cikojevic, L. V. Markic and J. Boronat, New J. Phys. 22, 053045 (2020).
  * (19) V. Cikojevic, L. V. Markic, M. Pi, M. Barranco and J. Boronat, Phys. Rev. A 102, 033335 (2020).
  * (20) I. Bloch, J. Dalibard and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * (21) T. Kinoshita, T. Wenger and D. S. Weiss, Science 305, 1125 (2004).
  * (22) B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, J. I. Cirac, G. V. Shlyapnikov, T. W. Hänsch and I. Bloch, Nature 429, 277 (2004).
  * (23) T. Giamarchi, Quantum Physics in One Dimension (Oxford Univ. Press, 2004).
  * (24) M. Kunitski, Q. Guan, H. Maschkiwitz, J. Hahnenbruch, S. Eckart, S. Zeller, A. Kalinin, M. Schöffler, L. Ph. H. Schmidt, T. Jahnke, D. Blume and R. Dörner, Nature Phys. (2020). https://doi.org/10.1038/s41567-020-01081-3.
  * (25) C. Chin, R. Grimm, P. Julienne and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * (26) M. Olshanii, Phys. Rev. Lett. 81, 938 (1998).
  * (27) P. F. Bedaque and U. van Kolck, Annu. Rev. Nucl. Part. Sci. 52, 339 (2002).
  * (28) H.- W. Hammer, A. Nogga and A. Schwenk, Rev. Mod. Phys. 85, 197 (2013).
  * (29) A. Bulgac, Phys. Rev. Lett. 89, 050402 (2002).
  * (30) M. Valiente, Phys. Rev. A 100, 013614 (2019).
  * (31) A. Pricoupenko and D. S. Petrov, e-print arXiv:2012.10429.
  * (32) R. Dashen, S. Ma and H. J. Bernstein, Phys. Rev. 187, 345 (1969).
  * (33) V. Pastukhov, Phys. Lett. A 383, 894 (2019).
  * (34) M. Valiente and V. Pastukhov, Phys. Rev. A 99, 053607 (2019).
  * (35) M. Valiente, Phys. Rev. A 103, L021302 (2020).
  * (36) M. Valiente, Phys. Rev. A 102, 053304 (2020).
  * (37) Y. Sekino and Y. Nishida, Phys. Rev. A 97, 011602 (2018).
  * (38) Y. Hou and J. E. Drut, Phys. Rev. Lett. 125, 050403 (2020).
  * (39) A. J. Czejdo, J. E. Drut, Y. Hou, J. R. McKenney and K. J. Morrell, Phys. Rev. A 101, 063630 (2020).
  * (40) Y. Hou, A. J. Czejdo, J. DeChant, C. R. Shill and J. E. Drut, Phys. Rev. A 100, 063627 (2019).
  * (41) J. E. Drut, J. R. McKenney, W. S. Daza, C. L. Lin and C. R. Ordóñez, Phys. Rev. Lett. 120, 243002 (2018).
  * (42) A. Pricoupenko and D. S. Petrov, Phys. Rev. A 100, 042707 (2019).
  * (43) G. Guijarro, A. Pricoupenko, G. E. Astrakharchik, J. Boronat and D. S. Petrov, Phys. Rev. A 97, 061605 (2018).
  * (44) Y. Nishida, Phys. Rev. A 97, 061603 (2018).
  * (45) L. Pricoupenko, Phys. Rev. A 97, 061604 (2018).
  * (46) L. Pricoupenko, Phys. Rev. A 99, 012711 (2019).
  * (47) S. R. Beane, Phys. Rev. A 82, 063610 (2010).
  * (48) S. R. Beane, Eur. Phys. J. D 72, 55 (2018).
  * (49) J. MacGuire, J. Math. Phys. 7, 123 (1966).
  * (50) F. Calogero and A. Degasperis, Phys. Rev. A 11, 265 (1975).
  * (51) B. Gertjerenken, T. P. Billam, L. Khaykovich, and C. Weiss , Phys. Rev. A 86, 033608 (2012).
  * (52) G. E. Astrakharchik, J. Boronat, I. L. Kurbakov, Yu. E. Lozovik and F. Mazzanti, Phys. Rev. A 81, 013612 (2010).
  * (53) M. Valiente and D. Petrosyan, J. Phys. B: At. Mol. Opt. Phys. 42, 121001 (2009).
  * (54) M. Lüscher, Comm. Math. Phys. 105, 153 (1986).
  * (55) J. R. Taylor, Scattering Theory. The Quantum Theory of Nonrelativistic Collisions (Dover, 2000).
  * (56) M. Lüscher, Comm. Math. Phys. 104, 177 (1986).
  * (57) S. König, D. Lee and H.- W. Hammer, Ann. Phys. 327, 1450 (2012).
  * (58) S. Stringari and J. Treiner, J. Chem. Phys. 87, 5021 (1987).
  * (59) R. A. Aziz, V. P. S. Nain, J. S. Carley, W. L. Taylor and G. T. McConville, J. Chem. Phys. 70, 4330 (1979).
  * (60) K. Polejaeva and A. Rusetsky, Eur. Phys. J. A 48, 67 (2012).
  * (61) K. Huang, Statistical Mechanics (Wiley, 1987).
  * (62) Y. Yan and D. Blume, Phys. Rev. Lett. 116, 230401 (2016).
  * (63) M. Valiente and N. T. Zinner, Few-body Syst. 56, 845 (2015).
  * (64) M. Valiente and D. Petrosyan, J. Phys. B: At. Mol. Opt. Phys. 41, 161002 (2008).
  * (65) S. König and D. Lee, Phys. Lett. B 779, 9 (2018).
  * (66) D. S. Petrov and G. E. Astrakharchik, Phys. Rev. Lett. 117, 100401 (2016).
  * (67) See supplemental material.

Supplemental Material: Quantum liquids and droplets with low-energy
interactions in one dimension

##  I Decoupling approximation

Here, we derive the approximate Hamiltonian (4) by using a decoupling
fluctuation expansion. This consists of approximating the product
$\delta(x_{i}-x_{j})\delta(x_{j}-x_{k})$ as

$\displaystyle\delta(x_{i}-x_{j})\delta(x_{j}-x_{k})$
$\displaystyle\approx\langle\delta(x_{i}-x_{j})\rangle\delta(x_{j}-x_{k})$
$\displaystyle+\langle\delta(x_{i}-x_{k})\rangle\delta(x_{i}-x_{j})$
$\displaystyle+\langle\delta(x_{j}-x_{k})\rangle\delta(x_{i}-x_{k})$
$\displaystyle-2\langle\delta(x_{i}-x_{j})\rangle\langle\delta(x_{j}-x_{k})\rangle.$
(S1)

For identical bosons, we have
$\Gamma\equiv\langle\delta(x_{i}-x_{j})\rangle=\langle\delta(x_{i^{\prime}}-x_{j^{\prime}})\rangle$
for all $i\neq j$ and $i^{\prime}\neq j^{\prime}$. Hamiltonian (1) reduces
immediately to Hamiltonian (4).

The ground state, for attractive interactions ($G_{N}(\Gamma)<0$), gives for
the energy $E_{0}(N)$ MacGuire1966S

$E_{0}(N)=-\frac{m\left[G_{N}(\Gamma)\right]^{2}}{24\hbar^{2}}N(N+1)(N-1)+C_{N}(\Gamma).$
(S2)

Using the Hellmann-Feynman theorem, and eliminating $\Gamma$, we obtain

$\displaystyle\Gamma$ $\displaystyle=-\frac{mG_{N}}{6\hbar^{2}}(N+1),$ (S3)
$\displaystyle G_{N}$
$\displaystyle=-|g_{0}|\left[1-\frac{1}{1+\frac{6\hbar^{2}}{mg_{3}}[(N+1)(N-2)]^{-1}}\right].$
(S4)

Relation (8) in the limit $N\to\infty$ is obtained by inserting Eqs. (S3) and
(S4) into Eq. (S2) and taking the limit $N\to\infty$.

##  II Universal asymptotics

Here, we derive Eq. (10) for the asymptotic behaviour of the density profile
of quantum droplets at large distances from the center. We begin by using the
results of Ref. Koenig2012S , particularized to one spatial dimension. When
one of the particles (e.g. particle 1) is asymptotically far from the rest of
the system, the ground state wave function $\psi_{N}$ factorizes as

$\displaystyle\psi_{N}$
$\displaystyle\to\gamma_{N}\sqrt{\kappa_{1,N}r_{1,N-1}}K_{-1/2}(\kappa_{1,N}r_{1,N-1})$
$\displaystyle\times\psi_{N-1}(x_{2},\ldots,x_{N}).$ (S5)

Above, $\gamma_{N}$ is the asymptotic normalization coefficient in the
$N\leftrightarrow(N-1)+1$ breakup channel,
$\hbar^{2}\kappa_{1,N}^{2}/2m=|E_{N}-E_{N-1}|(N-1)/N$ represents the binding
energy with respect to the ground state with one fewer particle, $K_{\nu}$ is
the modified Bessel function of the second kind of order $\nu$,
$r_{1,N-1}=|x_{1}-X_{N-1}|$, with $X_{N-1}=\sum_{i=2}^{N}x_{i}/(N-1)$ the
center of mass coordinate of the $(N-1)$-body subsystem, and $\psi_{N-1}$ is
its ground state wave function. The density profile $\rho_{N}(x)$ of the bound
state is given by Eq. (6). At long distance away from the bulk of the droplet,
we can introduce Eq. (S5) into Eq. (6), obtaining, as $|x|\to\infty$,

$\rho_{N}(x)\to\frac{\pi}{2}|\gamma_{1,N}|^{2}A_{N-1}\frac{N}{N-1}e^{-\frac{N}{N-1}2\kappa_{1,N}|x|},$
(S6)

where $A_{N-1}=\int\mathrm{d}^{N-2}z|\psi_{N-1}|^{2}$ is a model-independent
normalization factor. Collectively, $\mathbf{z}$ is an $(N-2)$-dimensional
vector containing all degrees of freedom of the $(N-1)$-particle system except
for the center-of-mass coordinate.

## References

  * [1] J. MacGuire, J. Math. Phys. 7, 123 (1966).
  * [2] S. König, D. Lee and H.- W. Hammer, Ann. Phys. 327, 1450 (2012).


# Two more ways of spelling Gini Coefficient with Applications

Marta Milewska  Centrum Wiskunde & Informatica (CWI), P.O. Box 94079, 1090 GB
Amsterdam, The Netherlands  Remco van der Hofstad  Department of Mathematics
and Computer Science, Eindhoven University of Technology, 5600 MB Eindhoven,
The Netherlands  Bert Zwart  Centrum Wiskunde & Informatica (CWI), P.O. Box
94079, 1090 GB Amsterdam, The Netherlands  Department of Mathematics and
Computer Science, Eindhoven University of Technology, 5600 MB Eindhoven, The
Netherlands

###### Abstract

We describe two new ways of defining and interpreting the Gini coefficient. We
investigate the potential application of the Gini coefficient to
epidemiological models focusing on the negative binomial distribution. To this
end, we derive several asymptotic results for the Gini coefficient. Finally,
applying both new representations, we compute the Gini coefficient for a few
probability distributions justifying the statement that the Gini coefficient
can be used also for models relying on distributions other than negative
binomial.

Keywords: Gini coefficient; excess distribution; epidemiology; superspreading
event; asymptotic expansions
AMS subject classification: 62G32, 92D30

##  1 Introduction

Variability is one of the most fundamental features in data analysis.
Statistics provides us with many tools suitable for measuring variability, the
most popular of which is the variance. However, it can be argued that in some
cases, for example when dealing with a distribution that significantly
deviates from the Gaussian curve, variance misses to convey some important
information. Therefore, different metric have been proposed, which also
provide quantitative information if the variance is infinite.

A particular metric known as the Gini mean difference (abbreviated to GMD),
first proposed in 1912 by the Italian statistician, Corrado Gini. Since then
GMD and the parameters that can be derived from it (such as the Gini
coefficient) have received quite some interest. However, they are mainly
applied to quantify and compare the income disparities among countries while
it is being suggested that this application does not reflect the entire
potential of the GMD.

Out of the Gini family, the most well-known index is the Gini coefficient,
again mostly used in the area of income inequalities. It was developed
independently of the GMD and defined with respect to the Lorenz curve
(explained in more detail further), which definition until today remains the
most well-known formulation of the Gini coefficient. In 1914, Corrado Gini
discovered and documented the relationship between the Gini coefficient and
the Gini mean difference which turns out to be very simple: the Gini
coefficient equals the GMD divided by twice the mean, assuming that the latter
is positive. As a result, the Gini coefficient ranges from 0 to 1, where 0
corresponds to perfect equality and 1 to the situation where only one
observation is positive.

So far at least 14 distinct alternative representations have been described,
each of them leading to a different interpretation. The readers interested in
investigating these various representations further should refer to [17] or
for a more concise guide to [20].

Our motivation for considering the Gini coefficient arises from the field of
epidemiology, where it has been realized that, apart from the mean
reproduction number, variability can also play an important role. Viruses such
as COVID-19 are considered to be highly overdispersed (see [5] and [7]). So
far, describing overdispersion was restricted to a particular parametric model
(see [10]). We want to make the case that the Gini coefficient provides a
measure of overdispersion that is not limited to any particular probability
distribution, and is consistent with the insights provided by existing
parametric models.

As computing the Gini coefficient with its classic definition might by arduous
for some probability distributions, we first add two new representations: a
probabilistic and an analytical one. We then apply the newly invented
methodology to compute the Gini coefficient for a few classic probability
distributions. We link our findings with the potential application to
measuring overdispersion by giving more attention to the negative binomial
distribution, which is often used in epidemiology to model the number of
secondary cases per infected. In particular, we examine its asymptotic
behavior.

This paper is organised as follows: in Section 2 we derive two new
representations of the Gini coefficient: a probabilistic and an analytical
one. In Section 3 we apply the newly invented methodology to compute the Gini
coefficient for a few classic probability distributions: the Poisson
distribution in Section 3.1, the exponential distribution in Section 3.2, the
geometric distribution in Section 3.3, the Pareto distribution in Section 3.4,
the uniform distribution in Section 3.5 and the negative binomial distribution
in Section 3.6. When possible, we refer to other examples of computing the
Gini coefficient for the mentioned distributions. In Section 4 we focus on
applications of the Gini coefficient for the negative binomial distribution,
which is very often used for epidemiological modeling. We explain why the Gini
coefficient should be considered as a measure of overdispersion in infectious
diseases. In Theorems 4.2 and 4.3, we prove the asymptotic behavior of the
Gini coefficient of the negative binomial distribution for small and large
parameters, respectively. Lastly, in Section 5 we close with a conclusion.

##  2 New representations of the Gini coefficient

The Gini coefficient is originally defined based on the Lorenz curve, which
plots the proportion of the total income of the population ($y$-axis) that is
cumulatively earned by the bottom $x\%$ of the population. Thus, the line at
$45$ degrees represents perfect equality of incomes. The Gini coefficient can
then be thought of as the ratio of the area that lies between the line of
equality and the Lorenz curve over the total area under the line of equality.
Before its link with the Gini mean difference was discovered, the coefficient
was actually called a ‘concentration ratio’. However, the graphical
representation is not always that useful. In order to derive alternative
expressions of the Gini index, we first recall its well-known [17, 20]
definition in terms of the mean difference.

######  Definition 2.1 (Gini coefficient).

The stochastic representation of the Gini coefficient of a random variable $X$
with finite mean $\mathbf{E}[X]$ is

$G(X)=\frac{\mathbf{E}[{|X_{1}-X_{2}|}]}{2\mathbf{E}[X]},$ (1)

where $X_{1}$ and $X_{2}$ are independent copies of $X$.

###  2.1 Representation in terms of a excess random variable.

Formula (1) is the most popular representation of the Gini coefficient.
However, it turns out that for many probability distributions it yields a
rather unpleasant form. Hence, we propose a more instructive representation in
terms of excess random variable:

######  Definition 2.2 (excess random variable).

Let $X$ be a non-negative random variable with $0<\mathbf{E}[X]<\infty$. We
define $X^{*}$, the excess version of $X$, as

$\displaystyle\mathbf{P}(X^{*}\geq
x)=\frac{1}{\mathbf{E}[X]}\int_{x}^{\infty}\mathbf{P}(X>s)\hskip
2.84544pt\text{d}s,\hskip 11.38092ptx\in\mathbb{R}_{\geq 0}$ (2)

######  Theorem 2.3 (Gini coefficient in terms of excess variable).

For $X\geq 0$ and $X^{*}$ independent,

$\displaystyle G(X)=\mathbf{P}(X^{*}\geq X),$ (3)

where $X^{*}$ is a excess random variable defined in (2).

###### Proof.

By the definition of the Gini coefficient stated in (1),

$\displaystyle G(X)$
$\displaystyle=\frac{1}{2\mathbf{E}[X]}\int_{0}^{\infty}\int_{0}^{\infty}|y-x|\hskip
2.84544pt\text{d}F_{X}(y)\hskip 2.84544pt\text{d}F_{X}(x)$
$\displaystyle=\frac{1}{\mathbf{E}[X]}\int_{0}^{\infty}\int_{x}^{\infty}(y-x)\hskip
2.84544pt\text{d}F_{X}(y)\hskip 2.84544pt\text{d}F_{X}(x).$ (4)

We write
$\int_{x}^{\infty}(y-x)\text{d}F_{X}(y)=\int_{x}^{\infty}\int_{x}^{y}\text{d}z\hskip
2.84544pt\text{d}F_{X}(y)$ and since the integrands are non-negative we can
apply Fubini’s theorem to change the order of integration in the latter and
rearrange (2.1) as

$\displaystyle G(X)$
$\displaystyle=\frac{1}{\mathbf{E}[X]}\int_{0}^{\infty}\int_{x}^{\infty}\int_{z}^{\infty}\text{d}F_{X}(y)\hskip
2.84544pt\text{d}z\hskip 2.84544pt\text{d}F_{X}(x)$
$\displaystyle=\int_{0}^{\infty}\int_{x}^{\infty}\frac{\mathbf{P}(X>z)}{\mathbf{E}[X]}\hskip
2.84544pt\text{d}z\hskip 2.84544pt\text{d}F_{X}(x)=\mathbf{P}(X^{*}\geq X).$

∎

Denote $\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$. Theorem 2.3 is valid both for
discrete and continuous random variables. However, for an
$\mathbb{N}_{0}$-valued random variable, $X^{*}$ in (3) can be replaced by a
discrete version of excess random variable, $X^{*}_{d}$, with the probability
distribution

$\displaystyle\mathbf{P}(X^{*}_{d}=x)=\frac{\mathbf{P}(X>x)}{\mathbf{E}[X]},\hskip
14.22636ptx=0,1,2,3\ldots$

######  Proposition 2.4.

Take a random variable $X$ taking values in $\mathbb{N}_{0}$. Then,
$X^{*}\stackrel{{\scriptstyle d}}{{=}}X^{*}_{d}+Y$ with $Y\sim Unif[0,1]$,
independent of $X_{d}^{*}$. Consequently, $G(X)=\mathbf{P}(X^{*}_{d}\geq X)$.

###### Proof.

Observe first that $X^{*}_{d}+Y$ has a density function $f_{X^{*}_{d}+Y}$ such
that

$\displaystyle
f_{X^{*}_{d}+Y}(z)=\frac{\mathbf{P}(X>x)}{\mathbf{E}[X]}=p_{x}\hskip
14.22636pt\text{for}\hskip 5.69046ptz\in[x,x+1).$

We now differentiate (2) to obtain the density of $X^{*}$:

$\displaystyle f_{X^{*}}(z)$
$\displaystyle=\frac{\text{d}}{\text{d}z}\mathbb{P}(X^{*}\leq
z)=\frac{\mathbf{P}(X>z)}{\mathbf{E}[X]}=p_{x}\hskip 8.5359pt\text{for}\hskip
2.84544ptz\in[x,x+1).$

Thus, $X^{*}\stackrel{{\scriptstyle d}}{{=}}X^{*}_{d}+Y$. Next, observe that,
for all $\mathbb{N}_{0}$-valued $X$,

$\displaystyle\mathbf{P}(X^{*}\geq
X)=\sum_{x=0}^{\infty}\mathbf{P}(X=x)\mathbf{P}(X^{*}\geq
x)=\sum_{x=0}^{\infty}\mathbf{P}(X=x)\mathbf{P}(X^{*}_{d}+Y\geq x).$ (5)

Since we only consider $x\in\mathbb{N}_{0}$, it is true that

$\displaystyle\mathbf{P}(X^{*}_{d}+Y\geq x)=\mathbf{P}(X^{*}_{d}\geq x)\cdot
1+\mathbf{P}(X^{*}_{d}=x-1)\cdot\mathbb{P}(Y=1)=\mathbf{P}(X^{*}_{d}\geq x).$

Substituting this into (5) completes the proof of the second statement. ∎

Since using $X^{*}_{d}$ to compute the Gini index simplifies some
calculations, throughout the paper we will use this definition when working
with $\mathbb{N}_{0}$-valued random variables.

###  2.2 Representation in terms of the Fourier Transform.

By noticing that the Gini coefficient is a convolution at $0$ of the
distribution of $X$ and the tail distribution of $-X^{*}$, one can arrive at a
different representation using Fourier theory. Since Fourier theory is
different in the discrete and continuous world, we make a distinction between
the two cases.

####  2.2.1 Discrete integer-valued case.

The following theorem provides a representation of the Gini coefficient of a
non-negative integer-valued random variable.

######  Theorem 2.5 (Fourier representation of Gini coefficient in discrete
case).

Let $X$ be a nonnegative integer-valued random variable with probability mass
function $p_{X}(j)=\mathbf{P}(X=j)$ and Fourier transform
$\hat{p}_{X}(\theta)=\sum_{j=0}^{\infty}\mathbf{P}(X=j)e^{i\theta j}$. Then

$\displaystyle
G(X)=\frac{1}{4\pi\mathbf{E}[X]}\int_{0}^{2\pi}\frac{[1-\hat{p}_{X}(\theta)-\mathbf{E}[X](e^{-i\theta}-1)]\cdot\hat{p}_{X}(-\theta)}{1-\cos{\theta}}\hskip
2.84544pt\text{d}\theta.$ (6)

######  Remark 2.6.

Since $\hat{p}_{X}(\theta)\hat{p}_{X}(-\theta)=||\hat{p}_{X}(\theta)||^{2}$,
(6) can be alternatively written as

$\displaystyle
G(X)=\frac{1}{4\pi\mathbf{E}[X]}\int_{0}^{2\pi}\frac{[1-\mathbf{E}[X](e^{-i\theta}-1)]\cdot\hat{p}_{X}(-\theta)-||\hat{p}_{X}(\theta)||^{2}}{1-\cos{\theta}}\hskip
2.84544pt\text{d}\theta.$ (7)

###### Proof.

Write

$\displaystyle G(X)$
$\displaystyle=\sum_{j=0}^{\infty}\mathbf{P}(X=j)\mathbf{P}(X^{*}_{d}\geq
j)=\sum_{j=0}^{\infty}\mathbf{P}(-X=-j)\mathbf{P}(X^{*}_{d}\geq j).$ (8)

Now, if we define the functions $f_{X}$ and $g_{X}$ as

$\displaystyle f_{X}(j)=\mathbf{P}(-X=j)\hskip 8.5359pt\text{and}\hskip
8.5359ptg_{X}(j)=\mathbf{P}(X^{*}_{d}\geq j),$

then we notice that (8) is a convolution of these functions, say $(f_{X}\star
g_{X})(n)$, at $n=0$, so $G(X)=(f_{X}\star g_{X})(0),$ which can be rearranged
with the help of Fourier theory. Applying the inverse Fourier transform and
the convolution theorem we obtain

$\displaystyle(f_{X}\star g_{X})(n)$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\hat{f_{X}}(\theta)\cdot\hat{g_{X}}(\theta)\cdot
e^{-in\theta}\hskip 2.84544pt\text{d}\theta,$ (9)

Substituting $n=0$ in (9), we arrive at the new representation of $G(X)$ as

$\displaystyle
G(X)=\frac{1}{2\pi}\int_{0}^{2\pi}\hat{f}_{X}(\theta)\cdot\hat{g}_{X}(\theta)\hskip
2.84544pt\text{d}\theta.$ (10)

It remains to compute the Fourier transforms $\hat{f}_{X}(\theta)$ and
$\hat{g}_{X}(\theta)$. It holds that

$\displaystyle\hat{f}_{X}(\theta)=\sum_{j=-\infty}^{0}e^{i\theta
j}\mathbf{P}(X=-j)=\hat{p}_{X}(-\theta).$ (11)

Further,

$\displaystyle\hat{g}_{X}(\theta)=\sum_{j=0}^{\infty}e^{i\theta
j}\sum_{m=j}^{\infty}\mathbf{P}(X^{*}_{d}=m)=\sum_{m=0}^{\infty}\mathbf{P}(X^{*}_{d}=m)\frac{e^{i(m+1)\theta}-1}{e^{i\theta}-1},$

where the last equality follows from changing the order of summation and
substituting
$\frac{e^{i(m+1)\theta}-1}{e^{i\theta}-1}=\sum_{j=0}^{m}e^{i\theta j}$. Now,
after extracting $\frac{1}{e^{i\theta}-1}$ and splitting the sum we arrive at

$\displaystyle\hat{g}_{X}(\theta)$
$\displaystyle=\frac{e^{i\theta}}{e^{i\theta}-1}\sum_{m=0}^{\infty}\mathbf{P}(X^{*}_{d}=m)e^{i\theta
m}-\frac{1}{e^{i\theta}-1}\sum_{m=0}^{\infty}\mathbf{P}(X^{*}_{d}=m)$
$\displaystyle=\frac{e^{i\theta}\hat{p}_{X^{*}_{d}}(\theta)-\hat{p}_{X^{*}_{d}}(0)}{e^{i\theta}-1}=\frac{e^{i\theta}\hat{p}_{X^{*}_{d}}(\theta)-1}{e^{i\theta}-1}.$
(12)

We now compute $\hat{p}_{X^{*}_{d}}(\theta)$ as

$\displaystyle\hat{p}_{X^{*}_{d}}(\theta)$
$\displaystyle=\sum_{m=0}^{\infty}e^{im\theta}\frac{\mathbf{P}(X>m)}{\mathbf{E}[X]}=\frac{1}{\mathbf{E}[X]}\sum_{m=0}^{\infty}e^{im\theta}\sum_{l>m}^{\infty}\mathbf{P}(X=l)$
$\displaystyle=\frac{1}{\mathbf{E}[X]}\sum_{l=1}^{\infty}\mathbf{P}(X=l)\frac{e^{i\theta
l}-1}{e^{i\theta}-1}=\frac{1}{\mathbf{E}[X]}\sum_{l=0}^{\infty}\mathbf{P}(X=l)\frac{e^{i\theta
l}-1}{e^{i\theta}-1},$

where we again changed the order of summation and computed
$\sum_{m=0}^{l-1}e^{im\theta}$ to obtain the last equality. After extracting
$\frac{1}{e^{i\theta}-1}$ and splitting the sum we have

$\displaystyle\hat{p}_{X^{*}_{d}}(\theta)$
$\displaystyle=\frac{1}{\mathbf{E}[X](e^{i\theta}-1)}\bigg{(}\sum_{l=0}^{\infty}\mathbf{P}(X=l)e^{i\theta
l}-\sum_{l=0}^{\infty}\mathbf{P}(X=l)\bigg{)}$
$\displaystyle=\frac{\hat{p}_{X}(\theta)-1}{\mathbf{E}[X](e^{i\theta}-1)}.$
(13)

Substituting (2.2.1) into (2.2.1) yields

$\displaystyle\hat{g}_{X}(\theta)$
$\displaystyle=\frac{e^{i\theta}\hat{p}_{X^{*}_{d}}(\theta)-1}{e^{i\theta}-1}=\frac{e^{i\theta}\frac{\hat{p}_{X}(\theta)-1}{\mathbf{E}[X](e^{i\theta}-1)}-1}{e^{i\theta}-1}$
$\displaystyle=\frac{e^{i\theta}}{(e^{i\theta}-1)^{2}}\frac{\hat{p}_{X}(\theta)-1}{\mathbf{E}[X]}-\frac{1}{e^{i\theta}-1}.$

Finally, we rewrite

$\frac{e^{i\theta}}{(e^{i\theta}-1)^{2}}=\frac{1}{(e^{i\theta}-1)(1-e^{-i\theta})}=-\frac{1}{||e^{i\theta}-1||^{2}},$

so that

$\displaystyle\hat{g}_{X}(\theta)=\frac{1-\hat{p}_{X}(\theta)-\mathbf{E}[X](e^{-i\theta}-1)}{||e^{i\theta}-1||^{2}\mathbf{E}[X]}.$
(14)

We derive the final formula from Theorem 2.5 by substituting (11) and (14)
into (10) and noticing that
$||e^{i\theta}-1||^{2}=\sin(\theta)^{2}+(\cos(\theta)-1)^{2}=2[1-\cos(\theta)].$
∎

####  2.2.2 Continuous case.

The following theorem provides an alternative representation of the Gini
coefficient of a continuous non-negative random variable:

######  Theorem 2.7 (Fourier representation of Gini coefficient in continuous
case).

Let $X$ be a non-negative random variable with probability density function
$p_{X}(x)$ and Fourier transform
$\hat{p}_{X}(\theta)=\int_{0}^{\infty}p_{X}(x)e^{i\theta x}\hskip
2.84544pt\text{d}x$. Then

$\displaystyle
G(X)=\frac{1}{2\pi\mathbf{E}[X]}\int_{\mathbf{R}}\frac{\hat{p}_{X}(\theta)(1-\hat{p}_{X}(-\theta)-i\theta\mathbf{E}[X])}{\theta^{2}}\hskip
2.84544pt\mathrm{d}\theta.$ (15)

###### Proof.

The main idea in this proof is exactly the same as it was in the discrete
case, though we need to define functions $f_{X}$ and $g_{X}$ in a slightly
different way. Consider a convolution of $f_{X}(x),g_{X}(x)$ with

$\displaystyle f_{X}(x)=p_{X}(x)\hskip 5.69046pt\text{and}\hskip
5.69046ptg_{X}(x)=\mathbf{P}(-X^{*}\leq x)=\mathbf{P}(X^{*}\geq-x),$

where $g_{X}$ is defined for negative values of $x$. Applying analogous
reasoning as in the proof of Theorem 2.5, we get

$\displaystyle G(X)=\frac{1}{2\pi}\int_{\mathbf{R}}\hat{f}_{X}(\theta)\hskip
2.84544pt\hat{g}_{X}(\theta)\hskip 2.84544pt\text{d}\theta.$ (16)

We again compute the corresponding Fourier transforms. We have

$\displaystyle\hat{g}_{X}(\theta)$
$\displaystyle=\int_{-\infty}^{0}\mathbf{P}(X^{*}\geq-x)e^{i\theta x}\hskip
2.84544pt\text{d}x=\int_{0}^{\infty}\mathbf{P}(X^{*}\geq y)e^{-i\theta
y}\hskip 2.84544pt\text{d}y$ $\displaystyle=\int_{0}^{\infty}e^{-i\theta
y}\int_{y}^{\infty}p_{X^{*}}(z)\hskip 2.84544pt\text{d}z\hskip
2.84544pt\text{d}y=\frac{1}{\mathbf{E}[X]}\int_{0}^{\infty}e^{-i\theta
y}\int_{y}^{\infty}\mathbf{P}(X>z)\hskip 2.84544pt\text{d}z\hskip
2.84544pt\text{d}y$
$\displaystyle=\frac{i}{\theta\mathbf{E}[X]}\int_{0}^{\infty}\mathbf{P}(X>z)(e^{-i\theta
z}-1)\hskip 2.84544pt\text{d}z,$

where the last equality follows after changing the order of integration and
computing $\int_{0}^{z}e^{-i\theta y}\hskip
2.84544pt\text{d}y=\frac{i}{\theta}(e^{-i\theta z}-1)$. We substitute
$\mathbf{P}(X>z)=\int_{z}^{\infty}f_{X}(y)\hskip 2.84544pt\text{d}y$ and
change the order of integration again to obtain

$\displaystyle\hat{g}_{X}(\theta)=\frac{i}{\theta\mathbf{E}[X]}\int_{0}^{\infty}f_{X}(y)\hskip
2.84544pt\int_{0}^{y}(e^{-i\theta z}-1)\hskip 2.84544pt\text{d}z\hskip
2.84544pt\text{d}y.$

Note that $\int_{0}^{y}(e^{-i\theta z}-1)\hskip
2.84544pt\text{d}z=\frac{i}{\theta}(e^{-i\theta y}-1)-y$. Hence,

$\displaystyle\hat{g}_{X}(\theta)$
$\displaystyle=\frac{-1}{\theta^{2}\mathbf{E}[X]}\int_{0}^{\infty}f_{X}(y)e^{-i\theta
y}\hskip
2.84544pt\text{d}y+\frac{1}{\theta^{2}\mathbf{E}[X]}\int_{0}^{\infty}f_{X}(y)\hskip
2.84544pt\text{d}y$
$\displaystyle\qquad-\frac{i}{\theta\mathbf{E}[X]}\int_{0}^{\infty}yf_{X}(y)\hskip
2.84544pt\text{d}y,$

which simplifies to

$\displaystyle\hat{g}_{X}(\theta)=\frac{1-\hat{p}_{X}(-\theta)-i\theta\mathbf{E}[X]}{\theta^{2}\mathbf{E}[X]}.$
(17)

We obtain (15) by substituting (17) in (16) noting that
$\hat{f}_{X}(\theta)=\hat{p}_{X}(\theta).$ ∎

##  3 Gini coefficient for various distributions

In this section, we apply our expressions to the computation of $G$ for
various distributions. In some cases, the examples are new, and sometimes, our
expressions provide easier derivations.

###  3.1 Poisson distribution

Take $X\sim\text{Poisson}(\lambda)$. The characteristic function equals
$\hat{p}_{X}(\theta)=\exp(\lambda(e^{i\theta}-1)).$ Consequently,
$\hat{p}_{X}(-\theta)=\exp(\lambda(e^{-i\theta}-1))$ and we substitute these
results into (6) to obtain

$\displaystyle
G(X)=\frac{1}{4\lambda\pi}\int_{0}^{2\pi}\frac{[1-\exp(\lambda(e^{i\theta}-1))-\lambda(e^{-i\theta}-1)]\cdot\exp(\lambda(e^{-i\theta}-1))}{1-\cos(\theta)}\hskip
2.84544pt\text{d}\theta.$

Alternatively, we could use (7). We compute

$\displaystyle||\hat{p}_{X}(\theta)||^{2}$
$\displaystyle=||\exp(\lambda(e^{i\theta}-1))||^{2}$
$\displaystyle=e^{-2\lambda}||e^{\lambda(\cos{\theta}+i\sin{\theta})}||^{2}=e^{-2\lambda(1-\cos{\theta})}.$

Substituting yields

$\displaystyle
G(X)=\frac{1}{4\lambda\pi}\int_{0}^{2\pi}\frac{(1-\lambda(e^{-i\theta}-1))e^{\lambda(e^{-i\theta}-1)}-e^{-2\lambda(1-\cos{\theta})}}{1-\cos(\theta)}\hskip
2.84544pt\text{d}\theta.$

The excess formula $G(X)=\mathbb{P}(X^{*}_{d}\geq X)$ does not seem to be
convenient in this example. An alternative representation of $G(X)$ in terms
of Bessel functions can be found in [16].

###  3.2 Exponential distribution

######  Proposition 3.1.

If $X\sim\text{Exp}(\lambda)$ then $G(X)=\frac{1}{2}$.

###### Proof.

Note that $X^{*}\stackrel{{\scriptstyle d}}{{=}}X$ and apply Theorem 2.3. ∎

###  3.3 Geometric distribution

Take $X\sim$ Geom$(p)$ with $\mathbf{P}(X=j)=(1-p)^{j}p$. It can be shown that
$\mathbf{P}(X^{*}_{d}\geq j)=(1-p)^{j}$. We plug it in $G(X)$ to arrive at

$\displaystyle
G(X)=p\sum_{j=0}^{\infty}(1-p)^{j}(1-p)^{j}=\frac{p}{p(2-p)}=\frac{1}{2-p}.$
(18)

In case of geometric distribution Fourier representation is also easy to
compute and yields the same result as (18).

Alternatively, one can consider a shifted version of the Geometric
distribution following [8]: take $Y\sim$ Geom’$(p)$ with
$\mathbf{P}(Y=j)=p(1-p)^{j-1}$ supported on $\\{1,2,3,...\\}$. It follows that
$\mathbf{P}(Y^{*}_{d}\geq j)=(1-p)^{j}$ and thus
$G(Y)=\frac{p}{1-p}\sum_{j=1}^{\infty}(1-p)^{2j}=\frac{1-p}{2-p}.$

###  3.4 Pareto distribution

Take $X\sim$ Pareto$(\alpha)$ with $f_{X}(x)=\frac{\alpha
x_{m}^{\alpha}}{x^{\alpha+1}}=$, $\mathbf{E}X=\frac{\alpha x_{m}}{\alpha-1}$
and assume $\alpha>1$. We first compute $\mathbf{P}(X^{*}\geq x)$ for $x\geq
x_{m}$:

$\displaystyle\mathbf{P}(X^{*}\geq x)$
$\displaystyle=\frac{1}{\mathbf{E}[X]}\int_{y=x}^{+\infty}\mathbf{P}(X>y)\hskip
2.84544pt\text{d}y=\frac{\alpha-1}{\alpha
x_{m}}\int_{y=x}^{+\infty}\bigg{(}\frac{x_{m}}{y}\bigg{)}^{\alpha}\hskip
2.84544pt\text{d}y$
$\displaystyle=\frac{\alpha-1}{\alpha}x_{m}^{\alpha-1}\bigg{[}\frac{x^{1-\alpha}}{1-\alpha}\bigg{]}_{x}^{+\infty}=\frac{1}{\alpha}\bigg{(}\frac{x_{m}}{x}\bigg{)}^{\alpha-1}.$

Note that $\log\big{(}\frac{X}{x_{m}}\big{)}\sim Exp(\alpha)$,
$\log\big{(}\frac{X^{*}}{x_{m}}\big{)}|X^{*}\geq x_{m}\sim Exp(\alpha-1)$.
Since $\mathbf{P}(Y\geq
Y^{\prime})=\frac{\beta^{\prime}}{\beta+\beta^{\prime}}$, where $Y\sim
Exp(\beta)$ and $Y^{\prime}\sim Exp(\beta^{\prime})$, we obtain

$\displaystyle\mathbf{P}(X^{*}\geq X)$ $\displaystyle=\mathbf{P}(X^{*}\geq
X|X^{*}\geq x_{m})\mathbf{P}(X^{*}\geq x_{m})$
$\displaystyle=\mathbf{P}(\log\bigg{(}\frac{X^{*}}{x_{m}}\bigg{)}\geq\log\bigg{(}\frac{X}{x_{m}}\bigg{)}|X^{*}\geq
x_{m})\mathbf{P}(X^{*}\geq x_{m})$
$\displaystyle=\frac{\alpha}{\alpha+\alpha-1}\cdot\frac{1}{\alpha}=\frac{1}{2\alpha-1}.$

This is not a new result, but its derivation is quite simple– compare e.g. [8,
Example 1].

###  3.5 Continuous uniform distribution

We have not been able to find this explicit example in the literature.

######  Proposition 3.2.

If $X\sim\mathcal{U}_{[a,b]}$, then $G(X)=\frac{b-a}{3(b+a)}$ =
$\frac{2}{b-a}\cdot\frac{{\mathrm{Var}}X}{\mathbf{E}X}$.

###### Proof.

$\displaystyle G(X)$ $\displaystyle=\int\mathbf{P}(X^{*}\geq x)f_{X}(x)\hskip
2.84544pt\text{d}x=\int_{a}^{b}\mathbf{P}(X^{*}\geq x)\frac{1}{b-a}\hskip
2.84544pt\text{d}x.$

Applying (2), substituting $\mathbf{E}X=\frac{a+b}{2}$ and changing the order
of integration yields

$\displaystyle G(X)$
$\displaystyle=\frac{2}{(a+b)(b-a)}\int_{a}^{\infty}(y-a)\mathbf{P}(X>y)\hskip
2.84544pt\text{d}y$
$\displaystyle=\frac{2}{(a+b)(b-a)^{2}}\int_{a}^{b}(y-a)(b-y)\hskip
2.84544pt\text{d}y$
$\displaystyle=\frac{2}{(a+b)(b-a)^{2}}\cdot\frac{1}{6}(b-a)^{3}=\frac{b-a}{3(b+a)}.$

∎

###  3.6 Negative binomial distribution

Let $X\sim NB(k,p)$ with probability mass function
$\mathbf{P}(X=j)=\binom{k+j-1}{j}p^{k}(1-p)^{j}$ and
$\mathbf{E}X=\frac{k(1-p)}{p}$. The classic definition of the Gini index (1)
yields

$\displaystyle
G(X)=\frac{1}{p}\sum_{i=0}^{\infty}(-1)^{i}\binom{k+i}{i}\Bigg{(}\frac{1-p}{p^{2}}\Bigg{)}^{i}\frac{(2i)!}{i!(i+1)!}.$
(19)

For a derivation see [16], formula (2.11) for the mean absolute difference of
the negative binomial; then divide by twice the mean to obtain (19). The
complexity of this expression was in fact our main motivation to develop new
representations. Because of the alternating sign the formula is not very
stable numerically and thus not handy to work with. The representation in
terms of the excess random variable does not seem very helpful either. The
computation leads to a rather non-intuitive formula containing three infinite
sums. Instead, we apply Theorem 2.5 to obtain the following representation:

######  Theorem 3.3.

Take $X\sim NB(k,p)$ with probability mass function
$\mathbf{P}(X=j)=\binom{k+j-1}{j}p^{k}(1-p)^{j}$ and
$\mathbf{E}X=\frac{k(1-p)}{p}$, where $p$ is the probability of success. Then

$\displaystyle G(X)=$
$\displaystyle\frac{p}{k(1-p)}\frac{1}{2\pi}\int_{0}^{2\pi}\frac{[1-\frac{k(1-p)}{p}(e^{-i\theta}-1)]\cdot\bigg{(}\frac{pe^{i\theta}}{e^{i\theta}+p-1}\bigg{)}^{k}-\bigg{(}\frac{p^{2}}{p^{2}+2(p-1)(\cos{\theta}-1)}\bigg{)}^{k}}{||e^{i\theta}-1||^{2}}\hskip
2.84544pt\text{d}\theta.$ (20)

######  Remark 3.4.

Note that since $\binom{k+j-1}{j}$ can be interpreted as
$\frac{\Gamma(k+j)}{\Gamma(j+1)\Gamma(k)}$, (3.3) also applies to non-integer
$k$.

###### Proof.

We want to apply Theorem 2.5 to the negative binomial variable $X$. Since

$\displaystyle\hat{p}_{X}(\theta)=p^{k}(1-(1-p)e^{i\theta})^{-k}=\frac{p^{k}}{(pe^{i\theta}-e^{i\theta}+1)^{k}},$
(21)

then

$\displaystyle\hat{p}_{X}(-\theta)=\frac{p^{k}}{(pe^{-i\theta}-e^{-i\theta}+1)^{k}}=\bigg{(}\frac{pe^{i\theta}}{e^{i\theta}+p-1}\bigg{)}^{k}.$
(22)

Substituting (21) and (22) into (6) yields

$\displaystyle G(X)$
$\displaystyle=\frac{p}{k(1-p)}\frac{1}{2\pi}\int_{0}^{2\pi}\frac{[1-\frac{p^{k}}{(pe^{i\theta}-e^{i\theta}+1)^{k}}-\frac{k(1-p)}{p}(e^{-i\theta}-1)]\cdot\bigg{(}\frac{pe^{i\theta}}{e^{i\theta}+p-1}\bigg{)}^{k}}{||e^{i\theta}-1||^{2}}\hskip
2.84544pt\text{d}\theta.$ (23)

We investigate the product of $\hat{p}_{X}(\theta)$ and
$\hat{p}_{X}(-\theta)$:

$\displaystyle\bigg{(}\frac{p}{pe^{i\theta}-e^{i\theta}+1}\cdot\frac{pe^{i\theta}}{e^{i\theta}+p-1}\bigg{)}^{k}=\bigg{(}\frac{p^{2}e^{i\theta}}{(pe^{i\theta}-e^{i\theta}+1)(e^{i\theta}+p-1)}\bigg{)}^{k}.$

However, we have that

$\displaystyle(pe^{i\theta}-$ $\displaystyle
e^{i\theta}+1)(e^{i\theta}+p-1)=pe^{2i\theta}+p^{2}e^{i\theta}-2pe^{i\theta}-e^{2i\theta}+2e^{i\theta}+p-1$
$\displaystyle=e^{i\theta}\big{(}p^{2}-2p+2+p(e^{i\theta}+e^{-i\theta})-(e^{i\theta}+e^{-i\theta})\big{)}$
$\displaystyle=e^{i\theta}\big{(}p^{2}+2(p-1)(\cos(\theta)-1)\big{)}.$

Therefore

$\displaystyle\frac{p^{2}e^{i\theta}}{(pe^{i\theta}-e^{i\theta}+1)(e^{i\theta}+p-1)}=\frac{p^{2}}{p^{2}+2(p-1)(\cos{\theta}-1)}.$

After plugging in the above computation in (3.6) and rearranging we arrive at
the expression (3.3). ∎

In the next section, we discuss the asymptotic behavior of the Gini
coefficient of the negative binomial distribution as $k\to 0$ and
$k\to\infty$.

##  4 Applications to epidemics

In epidemiology, the basic reproduction number $R_{0}$ denotes the expected
number of new infections directly generated by one case. If $R_{0}>1$, the
epidemic is growing, and the epidemic is dying our when $R_{0}<1$. It has been
determined that transmission patterns of SARS-COV-2 and many other viruses are
highly overdispersed (see [5] or [7]). Looking at $R_{0}$ alone might
therefore be misleading. In case of COVID-19 in particular, it is estimated
that more than 80$\%$ of new infections resulted from the top 10$\%$ of most
infectious individuals (see [1], [14] or [15]). This means that the vast
majority of infected individuals will not pass the infection on to anyone.
$R_{0}$ should therefore be complemented with another statistic that sheds
light on possible overdispersion in the number of secondary cases.

To this end, the authors in [10] assume a parametric setting, in which the
number of secondary cases follows a negative binomial distribution with
parameters $k$ and $R_{0}$, which results in its variance being equal to
$R_{0}(1+R_{0}/k)$. In such a setting, the smaller the value of $k$, the
greater heterogeneity in the distribution of secondary infections. Despite
coming from a very particular model, this interpretation became so popular
that multiple papers adopted the symbol $k$ and referred to it as an
‘overdispersion parameter’ (see for instance [1], [4], [9] or [11]). Even non-
mathematical press followed this trend, calling $k$ the key to overcoming the
pandemic of COVID-19 (see [18]).

In this section, we investigate the relation between $k$ and the Gini index.
We show numerically that the Gini index, when specialized to the negative
binomial distribution, depends monotonically on $k$, and complement this with
rigorous asymptotic estimates for small and large $k$. This indicates that the
Gini index provides insights consistent with the ones obtained so far, without
the need to rely on a parametric model. Thus, the Gini index can serve as
candidate to measure the variability in the infectiousness of a disease.

This section is organised as follows: in Section 4.1 we first enable a
continuous interpretation of parameter $k$ of the negative binomial
distribution by expressing is as a composition of Poisson and Gamma processes.
Next we investigate monotonicity of the Gini coefficient in $k$. To that end,
we illustrate its behavior in plots and support them with statements about
asymptotics of the Gini coefficient for $k\to 0$ and $k\to\infty$. In Section
4.2 we show for the real-world data set that the Gini coefficient provides
insights consistent with using a parametric model based on $k$.

###  4.1 Asymptotic behavior of the Gini coefficient

In Section 3.6 we showed that the Gini index for the negative binomial
distribution cannot be easily obtained from (1), neither by its excess
representation $\mathbb{P}(X^{*}_{d}\geq X).$ In this section, we investigate
the Gini index for small and large values of $k$. For that purpose it is
useful to first rewrite the negative binomial distribution in terms of the
composition of Poisson and Gamma processes, which has the added benefit that
the resulting interpretation is natural for non-integer values of $k$.

Let $X$ follow a $Poisson(\nu)$ distribution, where the rate $\nu$ is also a
random variable and it follows a $Gamma(k,\lambda)$ distribution. To maintain
a clear connection to epidemiological notation, we denote
$\mathbf{E}[\nu]=R_{0}$.

######  Lemma 4.1.

For $X$ as above, it holds that $X\sim NB(k,p)$ with parameter $k$,
corresponding to the number of ‘successes’, and parameter
$p=\frac{\lambda}{\lambda+1}$ denoting the probability of success.

The proof of Lemma 4.1 is straightforward and is hence omitted. Note that
$p=(1+R_{0}/k)^{-1}$. In the following, we will write $\Gamma_{k}$ referring
to a Gamma process determined by $\Gamma(k,\lambda)$ and $N(\Gamma_{k})$ to
denote a subordinated Poisson process with rate following this Gamma process.

With the new, continuous interpretation in mind we will have a closer look at
the behaviour of the Gini coefficient as $k$ increases. For creating some
visual intuition, Figure 1 plots the Gini coefficient as a function of $k$ for
different values of $p=\frac{\lambda}{\lambda+1}$.

(a) $p=0.9$

(b) $p=0.5$

(c) $p=0.1$

Figure 1: Plots of the Gini coefficient of negative binomial distribution

Figure 1 suggests that the Gini coefficient of the negative binomial is
decreasing in $k$. Unfortunately, a rigorous proof of the statement remains an
open problem. However, we succeeded in showing that the limiting behaviour of
the Gini coefficient of the negative binomial distribution aligns with the
desired statement, i.e., $G(NB(k,p))$ approaches $1$ as $k$ approaches $0$ and
approaches $0$ as $k$ approaches infinity. We formulate these results in the
following two theorems:

######  Theorem 4.2 (Asymptotics of the Gini coefficient for $k$ small).

Let $X\sim NB(k,p)$. Then, for $k\to 0$,

$\displaystyle G(NB(k,p))=1+c\cdot k+o(k),$ (24)

with
$c=2\log(p)+\frac{2-p}{1-p}\log\big{(}2p(1-p)\big{)}=2\log\bigg{(}\frac{\lambda}{\lambda+1}\bigg{)}-(\lambda+2)\log\bigg{(}\frac{\lambda^{2}+2\lambda}{(\lambda+1)^{2}}\bigg{)}$.

######  Theorem 4.3 (Asymptotic behavior of the Gini coefficient for $k$
large).

Let $X\sim NB(k,p)$. Then, for $k\to\infty$,

$\displaystyle\sqrt{k}\hskip
2.84544ptG(NB(k,p))\longrightarrow\frac{\sqrt{1+\lambda}}{\sqrt{\pi}}.$ (25)

(a) $p=0.9$

(b) $p=0.5$

(c) $p=0.1$

Figure 2: Behavior of the Gini coefficient for $k\to 0$ vs. $y=1+ck$

(a) $p=0.9$

(b) $p=0.5$

(c) $p=0.1$

Figure 3: Behavior of the Gini coefficient times $\sqrt{k}$ for $k\to\infty$
vs. $y=\frac{\sqrt{1+\lambda}}{\sqrt{\pi}}$

Proofs of these theorems are given in the Appendix. We validate the asymptotic
estimates in the theorems by plotting the Gini coefficient and its limiting
functions against each other. See Figure 2 for $k\to 0$ and Figure 3 for
$k\to\infty$.

###  4.2 The Gini coefficient for COVID-19’ secondary infections

In this section we investigate how our formulas perform in a real-world
example. Heterogeneity in the number of secondary cases receives significant
attention during the currently ongoing pandemic of COVID-19. [11] aims to
quantify overdispersion by estimating $k$. We use the data provided there to
first compute the Gini coefficient explicitly (by (1)) and then calculate it
using formulas we derived in this paper. For the latter, we substituted the
parameters $p$ and $k$ obtained by the authors. The table below presents all
results.

Location | Jakarta - Depok | Batan
---|---|---
$R_{0}$ | 6,79 | 2,47
$k$ | 0,06 | 0,2
$p$ | 0,008 | 0,06
$G(X)$ computed explicitly from data | 0,9213411 | 0,83191721
$G(X)$ computed with (6) | 0,9269144 | 0,8151066
$G(X)$ computed with (24) | 0.9193061 | 0.7623808
Table 1: Secondary cases data summary for Jakarta-Depok and Batan.

These results show that our formulas come very close to the actual value of
the Gini coefficient. It is also apparent that they reflect the dynamic of
epidemic - lower values of $k$ assigned to higher overdispersion correspond to
higher value of the Gini coefficient.

##  5 Conclusion

The models based on the negative binomial distribution of secondary infections
have been used to estimate overdispersion in the spread of infectious diseases
for a long time. However, this approach necessitates parametric models
suitable only for specific situations. We argue that it is preferable to use
the Gini coefficient to measure overdispersion. It correctly captures the
dynamic of epidemic and is arguably more informative than popular so far
parameter $k$ as it is solely known as the measure of variability and it is
not bound to any specific probability distribution.

Thanks to the new representations that we presented in this paper, we can
express the Gini coefficient as a simple function of the parameters of the
distribution for multiple probability distributions, as long as they have
finite expected value, making the Gini index a versatile statistic not
depending on parametric assumptions. In particular, it can be useful in the
case of heavy-tail distributions where the variance might be infinite. In [6]
the authors present methods that allow estimating the Gini coefficient for
heavy-tailed data with infinite variance.

The applications of the newly developed representations do not need to be
limited to epidemiology. The Gini coefficient was advocated as an indicator in
many other areas and we hope our results could help widen its usage. For
reference see for example [3], where the Gini coefficient is suggested as a
measure of heart rate variability to assess the level of mental stress, or
[19], where the Gini coefficient is applied to explain inequalities in
resource use around the globe.

##  Appendix A Asymptotics Gini coefficient of negative binomial

###  A.1 Proof of Theorem 4.2

We begin by proving two auxiliary propositions.

######  Proposition A.1 (Asymptotic behavior of $\mathbf{P}(X_{k}=j)$).

When $k\to 0$,

$\displaystyle\mathbf{P}(X_{k}=0)=1-k\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(k),$
(26) $\displaystyle\mathbf{P}(X_{k}=j)=c_{j}k+o(k),\hskip 28.45274ptj\geq 1,$
(27)

where $c_{j}=\frac{1}{j(\lambda+1)^{j}}.$ It is also true that
$\mathbf{E}[X_{k}]$ is asymptotically equal to $k$ for $k$ small.

###### Proof.

For $X_{k}\sim Poisson(\Gamma(k,\lambda))$,

$\displaystyle\mathbf{P}(X_{k}=j)$
$\displaystyle=\frac{\lambda^{k}}{j!\Gamma(k)}\frac{\Gamma(j+k)}{(\lambda+1)^{j+k}}=\bigg{(}\frac{\lambda}{\lambda+1}\bigg{)}^{k}\frac{1}{j!(\lambda+1)^{j}}\frac{\Gamma(j+k)}{\Gamma(k)}.$

Hence, for $k\to 0$,

$\displaystyle\mathbf{P}(X_{k}=0)$
$\displaystyle=\bigg{(}\frac{\lambda}{\lambda+1}\bigg{)}^{k}=e^{-k\log(\frac{\lambda+1}{\lambda})}=1-k\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(k),$

where we apply $e^{-ak}=1-ak+o(k)$ in the second equality. Further, applying
the same asymptotic formula together with $\Gamma(k)\sim\frac{1}{k}$ and
$\Gamma(j+k)=\Gamma(j)+O(k)$ (see [2, 6.1.35]), we obtain, for $j\geq 1$,

$\displaystyle\mathbf{P}(X_{k}=j)$
$\displaystyle=\bigg{(}1-k\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(k)\bigg{)}\frac{1}{j!(\lambda+1)^{j}}\big{(}k\Gamma(j)+o(k)\big{)}$
$\displaystyle=k\frac{(j-1)!}{j!(\lambda+1)^{j}}+o(k),$

which yields (27). ∎

######  Proposition A.2 (Asymptotic behavior of $\mathbf{P}(X_{k}^{*}\geq
j)$).

For $k\to 0$ and $j\geq 1$,

$\displaystyle\mathbf{P}(X_{k}^{*}\geq j)$ $\displaystyle=d^{*}_{j}+o(1),$
(28)

with
$d^{*}_{j}=\big{(}\frac{1}{\lambda+1}\big{)}^{j}-j\lambda\sum_{n=j}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}$.

###### Proof.

We first compute for $l\geq 1$ and $k\to 0$,

$\displaystyle\mathbf{P}(X^{*}_{k}=l)$
$\displaystyle=\frac{\mathbf{P}(X_{k}>l)}{\mathbf{E}X_{k}}=\frac{\lambda}{k}\mathbf{P}(X_{k}\geq
l+1)=\frac{\lambda}{k}\left(1-\sum_{m=0}^{l}\mathbf{P}(X_{k}=m)\right)$
$\displaystyle=$
$\displaystyle\frac{\lambda}{k}\bigg{(}1-\sum_{m=0}^{l}[c_{m}k+o(k)]\bigg{)}=\frac{\lambda}{k}\bigg{(}k\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(k)-\sum_{m=1}^{l}c_{m}k+o(k)\bigg{)}$
$\displaystyle=$
$\displaystyle\lambda\sum_{m=l+1}^{\infty}c_{m}+o(1)=d_{l}+o(1),$

where

$d_{l}=\lambda\sum_{m=l+1}^{\infty}\frac{1}{m(\lambda+1)^{m}}=\frac{\lambda}{(\lambda+1)^{l+1}}\sum_{n=0}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n}}{n+l+1},$

and in the fourth equality, we have substituted (27). Hence, in a similar
manner,

$\mathbf{P}(X^{*}_{k}\geq
j)=1-\mathbf{P}(X^{*}_{k}<j)=1-\sum_{l=0}^{j-1}d_{l}+o(1),$

and we compute

$\displaystyle\sum_{l=0}^{j-1}d_{l}$
$\displaystyle=\sum_{l=0}^{j-1}\frac{\lambda}{(\lambda+1)^{l+1}}\sum_{n=0}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n}}{n+l+1}=\lambda\sum_{l=0}^{j-1}\sum_{n=0}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+l+1}}{n+l+1}$
$\displaystyle=\lambda\sum_{l=0}^{j-1}\sum_{n=l}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}=j\lambda\sum_{n=j}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}+\lambda\sum_{n=0}^{j-1}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}(n+1)$
$\displaystyle=j\lambda\sum_{n=j}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}+\lambda\cdot\frac{\lambda+1-\big{(}\frac{1}{\lambda+1}\big{)}^{j-1}}{\lambda(\lambda+1)}$
$\displaystyle=1-\bigg{(}\frac{1}{\lambda+1}\bigg{)}^{j}+j\lambda\sum_{n=j}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1},$

from which the statement follows. ∎

Now we can proceed with the proof of the main theorem:

###### Proof of Theorem 4.2.

Since $\mathbf{P}(X^{*}_{k}\geq 0)=1$ we write, for all $N$,

$\displaystyle G(X_{k})$
$\displaystyle=\mathbf{P}(X_{k}=0)+\sum_{j=1}^{\infty}\mathbf{P}(X_{k}=j)\mathbf{P}(X^{*}_{k}\geq
j)$
$\displaystyle=1-k\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(k)+\sum_{j=1}^{N}(c_{j}k+o(k))(d^{*}_{j}+o(1))$
$\displaystyle+\sum_{j=N+1}^{\infty}\mathbf{P}(X_{k}=j)\mathbf{P}(X^{*}_{k}\geq
j),$

where we have substituted (26) and (27) from Proposition A.1 and (28) from
Proposition A.2. Thus, for all $N$

$\displaystyle\limsup_{k\to
0}\frac{1-G(X_{k})}{k}\leq\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(1)-\sum_{j=1}^{N}c_{j}d^{*}_{j}.$

and therefore, taking $N\to\infty$

$\displaystyle\limsup_{k\to
0}\frac{1-G(X_{k})}{k}\leq\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(1)-\sum_{j=1}^{\infty}c_{j}d^{*}_{j}.$
(29)

On the other hand, for all $N$ we have

$\displaystyle\liminf_{k\to 0}\frac{1-G(X_{k})}{k}$
$\displaystyle\geq\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(1)-\sum_{j=1}^{N}c_{j}d^{*}_{j}-\frac{1}{k}\sum_{N+1}^{\infty}\mathbf{P}(X_{k}=j)$
$\displaystyle\geq\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(1)-\sum_{j=1}^{N}c_{j}d^{*}_{j}-\frac{1}{k}\frac{\mathbf{E}[X_{k}]}{N+1}$
$\displaystyle=\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(1)-\sum_{j=1}^{N}c_{j}d^{*}_{j}-\frac{1}{\lambda(N+1)},$

where we have applied $\mathbf{P}(X_{k}=j)\leq 1$ and Markov inequality.
Hence, letting $N\to\infty$,

$\displaystyle\liminf_{k\to
0}\frac{1-G(X_{k})}{k}\geq\log\big{(}\frac{\lambda+1}{\lambda}\big{)}+o(1)-\sum_{j=1}^{\infty}c_{j}d^{*}_{j}.$
(30)

Combining (30) and (29) we can write

$\displaystyle
G(X_{k})=1+k\bigg{(}\sum_{j=1}^{\infty}c_{j}d^{*}_{j}-\log\big{(}\frac{\lambda+1}{\lambda}\big{)}\bigg{)}+o(k).$
(31)

Now we compute the term $\sum_{j=1}^{\infty}c_{j}d^{*}_{j}$. We have

$\displaystyle\sum_{j=1}^{\infty}c_{j}d^{*}_{j}$
$\displaystyle=\sum_{j=1}^{\infty}\frac{1}{j(\lambda+1)^{j}}\bigg{(}\bigg{(}\frac{1}{\lambda+1}\bigg{)}^{j}-j\lambda\sum_{n=j}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}\bigg{)}$
$\displaystyle=\sum_{j=1}^{\infty}\frac{1}{j(\lambda+1)^{2j}}-\lambda\sum_{j=1}^{\infty}\frac{1}{(\lambda+1)^{j}}\sum_{n=j}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}.$
(32)

The first summand in (A.1) is easily computable and equals
$-\log\big{(}\frac{\lambda^{2}+2\lambda}{(\lambda+1)^{2}}\big{)}$. The second
summand can be computed by interchanging the order of summation, as

$\displaystyle\lambda$
$\displaystyle\sum_{j=1}^{\infty}\frac{1}{(\lambda+1)^{j}}\sum_{n=j}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}=\lambda\sum_{n=1}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}\sum_{j=1}^{n}\frac{1}{(\lambda+1)^{j}}$
$\displaystyle=\sum_{n=1}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}\bigg{(}1-\big{(}\frac{1}{\lambda+1}\big{)}^{n}\bigg{)}=\sum_{n=1}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{n+1}}{n+1}-\sum_{n=1}^{\infty}\frac{\big{(}\frac{1}{\lambda+1}\big{)}^{2n+1}}{n+1}$
$\displaystyle=-\log\big{(}\frac{\lambda}{\lambda+1}\big{)}-\frac{1}{\lambda+1}+\big{(}\lambda+1\big{)}\log\bigg{(}\frac{\lambda^{2}+2\lambda}{(\lambda+1)^{2}}\bigg{)}+\frac{1}{\lambda+1}$
$\displaystyle=\big{(}\lambda+1\big{)}\log\bigg{(}\frac{\lambda^{2}+2\lambda}{(\lambda+1)^{2}}\bigg{)}-\log\big{(}\frac{\lambda}{\lambda+1}\big{)}.$

Hence, (31) becomes

$\displaystyle
G(X_{k})=1+k\bigg{(}2\log\bigg{(}\frac{\lambda}{\lambda+1}\bigg{)}-(\lambda+2)\log\bigg{(}\frac{\lambda^{2}+2\lambda}{(\lambda+1)^{2}}\bigg{)}\bigg{)}+o(k),$

what can be also written as

$\displaystyle
G(X_{k})=1+k\big{[}2\log(p)+\frac{2-p}{1-p}\log\big{(}2p(1-p)\big{)}\big{]}+o(k),$

when we take $p=\frac{\lambda}{\lambda+1}$. ∎

###  A.2 Proof of Theorem 4.3

We first prove the following proposition:

######  Proposition A.3.

Denote $X_{k}=N(\Gamma_{k})$. For $k\to\infty$,

$\displaystyle\frac{X_{k}-\frac{k}{\lambda}}{\sigma\sqrt{k}}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N(0,1)\hskip 14.22636pt\text{as}\hskip
8.5359ptk\to\infty,$

where $\sigma=\frac{\sqrt{1+\lambda}}{\lambda}$.

###### Proof.

We can write

$\displaystyle\frac{N(\Gamma_{k})-\frac{k}{\lambda}}{\sqrt{k}}=\frac{N(\Gamma_{k})-\Gamma_{k}}{\sqrt{\Gamma_{k}}}\cdot\sqrt{\frac{\Gamma_{k}}{k}}+\frac{\Gamma_{k}-\frac{k}{\lambda}}{\sqrt{k}}.$

Denote $\tilde{\Gamma}_{k}=\frac{\lambda\Gamma_{k}-k}{\sqrt{k}}$. From the
central limit theorem, it follows that

$\displaystyle\mathbf{P}(\tilde{\Gamma}_{k}\geq x)\to\bar{\Phi}(x)\hskip
14.22636pt\text{as}\hskip 8.5359ptk\to\infty.$

Hence,

$\displaystyle\frac{\Gamma_{k}-\frac{k}{\lambda}}{\sqrt{k}}=\frac{\lambda\Gamma_{k}-k}{\lambda\sqrt{k}}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}W_{\frac{1}{\lambda^{2}}}\hskip
14.22636pt\text{as}\hskip 8.5359ptk\to\infty.$

Furthermore, by the strong law of large numbers

$\displaystyle\frac{\Gamma_{k}}{k}\to\frac{1}{\lambda}\hskip
14.22636pt\text{as}\hskip 8.5359ptk\to\infty,$

and thus $\sqrt{\frac{\Gamma_{k}}{k}}\to\frac{1}{\sqrt{\lambda}}.$ Finally,
$\Gamma_{k}\to\infty$ a.s. as $k\to\infty$ and again, by the central limit
theorem we know that for a Poisson process $N(k)$ it is true that

$\displaystyle\frac{N(k)-k}{\sqrt{k}}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}W_{1}\hskip 14.22636pt\text{as}\hskip
8.5359ptk\to\infty.$

We conclude that also
$\frac{N(\Gamma_{k})-\Gamma_{k}}{\sqrt{\Gamma_{k}}}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}W_{1}$ as $k\to\infty$. The claim follows. ∎

###### Proof of Theorem 4.3.

Recall $X_{k}=N(\Gamma_{k})$ and set
$\tilde{X}_{k}=\frac{X_{k}-\frac{k}{\lambda}}{\sigma\sqrt{k}}$. From (1) we
obtain

$\displaystyle\sqrt{k}G(X_{k})=\frac{\mathbf{E}\big{[}{|X_{k}-X^{\prime}_{k}|}\big{]}}{2\mathbf{E}[X]}=\frac{\lambda}{2\sqrt{k}}\mathbf{E}\big{[}{|X_{k}-X^{\prime}_{k}|}\big{]},$

where $X^{\prime}_{k}$ is an independent copy of $X_{k}$. Then, for fixed $M$,

$\displaystyle\sqrt{k}G(X_{k})$
$\displaystyle=\frac{\lambda}{2\sqrt{k}}\mathbf{E}\big{[}{|X_{k}-X^{\prime}_{k}|}\cdot\mathbf{1}_{\\{|X_{k}-X^{\prime}_{k}|<M\sqrt{k}\\}}\big{]}$
$\displaystyle+\frac{\lambda}{2\sqrt{k}}\mathbf{E}\big{[}{|X_{k}-X^{\prime}_{k}|}\cdot\mathbf{1}_{\\{|X_{k}-X^{\prime}_{k}|>M\sqrt{k}\\}}\big{]}=I+II.$

We study parts $I$ and $II$ separately, starting with $I$:

$\displaystyle I$
$\displaystyle=\frac{\lambda\sigma\sqrt{k}}{2\sqrt{k}}\mathbf{E}\big{[}{|\tilde{X}_{k}-\tilde{X^{\prime}}_{k}|}\cdot\mathbf{1}_{\\{|\tilde{X}_{k}-\tilde{X^{\prime}}_{k}|<\frac{M\sqrt{k}}{\sigma\sqrt{k}}\\}}\big{]}$
$\displaystyle=\frac{\sigma\lambda}{2}\mathbf{E}\big{[}{|\tilde{X}_{k}-\tilde{X^{\prime}}_{k}|}\cdot\mathbf{1}_{\\{|\tilde{X}_{k}-\tilde{X^{\prime}}_{k}|<\frac{M}{\sigma}\\}}\big{]}$

Letting first $k$ and then $M$ tend to infinity, by Proposition A.3, we arrive
at

$\displaystyle
I\to\frac{\sigma\lambda}{2}\mathbf{E}\big{[}{|\tilde{X}-\tilde{X^{\prime}}|}\cdot\mathbf{1}_{\\{|\tilde{X}-\tilde{X^{\prime}}|<\infty\\}}\big{]},$

where $\tilde{X}$ and $\tilde{X^{\prime}}$ are independent standard normal
random variables. Write $Y=\tilde{X}-\tilde{X^{\prime}}$. Then $Y\sim N(0,2)$
and $\mathbb{E}[|Y|]=\frac{2}{\sqrt{\pi}}$. Thus,

$\displaystyle
I\to\frac{\sigma\lambda}{\sqrt{\pi}}=\frac{\sqrt{1+\lambda}}{\sqrt{\pi}}.$
(33)

For term $II$,

$\displaystyle II$
$\displaystyle=\frac{\lambda}{2\sqrt{k}}\mathbf{E}\big{[}{|X_{k}-X^{\prime}_{k}|}\cdot\mathbf{1}_{\\{|X_{k}-X^{\prime}_{k}|>M\sqrt{k}\\}}\big{]}\leq\frac{\lambda}{2\sqrt{k}}\frac{\mathbf{E}\big{[}{(X_{k}-X^{\prime}_{k})^{2}}\big{]}}{M\sqrt{k}}$
$\displaystyle=\frac{\lambda}{2\sqrt{k}}\frac{2{\mathrm{Var}}(X_{k})}{M\sqrt{k}}=\frac{\lambda}{2\sqrt{k}}\frac{2}{M\sqrt{k}}\frac{k(1+\lambda)}{\lambda^{2}}=\frac{1+\lambda}{M\lambda}\stackrel{{\scriptstyle
M\to\infty}}{{\longrightarrow}}0.$ (34)

Combining (33) and (A.2) we obtain the desired claim (25) from Theorem 4.3. ∎

##### Acknowledgement.

The work of MM is supported by the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement no.
945045, and by the NWO Gravitation project NETWORKS under grant no.
024.002.003.
The work of RvdH is supported in parts by the NWO through the Gravitation
Networks grant 024.002.003.

## References

  * [1] Abbot, S. Endo, A. et. al. (2020). Estimating the overdispersion in COVID-19 transmission using outbreak sizes outside China. _Wellcome Open Research_ 5: 67\.
  * [2] Abramowitz, M. and Stegun, I. (1964). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.
  * [3] Acosta-Batista C., Andreu-Heredia A., Carrazana-Escalona R., Cortina-Reyna S., Ghiya S., Saá-Muñoz N.A., Sánchez-Hechavarría M.E., (2019). Introduction of Application of Gini Coefficient to Heart Rate Variability Spectrum for Mental Stress Evaluation. _Arquivos Brasileiros de Cardiologia_ 113(4)
  * [4] Chen X., Chong M., Guo z., Huang Z., Wang J., Wang M., Wong E., yeoh E., Zee B., Zhao S. and Zhuang Z. (2021). Superspreading and heterogeneity in transmission of SARS, MERS, and COVID-19: A systematic review. _Computational and Structural Biotechnology Journal_ 19 5039–5046.
  * [5] Cirillo, P. and Taleb, N.N. (2020). Tail Risk of Contagious Diseases. _Nature Physics_ 16(6), 606–613.
  * [6] Cirillo, P. Fontanari, A. and Taleb, N.N. (2018). Gini estimation under infinite variance. _Physica A: Statistical Mechanics and its Applications_ 502, 256–269.
  * [7] Collins, J.J. and Wong, F. (2020). Evidence that coronavirus superspreading is fat-tailed. _Proceedings of the National Academy of Sciences_ 117(47).
  * [8] Dorfman, R. (1979). A Formula for the Gini Coefficient. _The Review of Economics and Statistics_ 61(1), 146–149
  * [9] Eilersen A., Nielsen B.F., Simonsen L. and Kim Sneppen K. (2021). Lockdowns exert selection pressure on overdispersion of SARS-CoV-2 variants. _medRxiv_
  * [10] Getz, W. M. Kopp, P. E. Lloyd-Smith, J. O. and Schreiber, S. J. (2005). Superspreading and the effect of individual variation on disease emergence. _Nature_ 438, 355–359.
  * [11] Hasan, A., Lestari, B., Kasim, M. F., Nuraini, N., Susanto, H. Triany, D. and Widyastuti, W. (2020). Superspreading in early transmissions of COVID‑19 in Indonesia. _Scientific Reports_ 10, 1–4.
  * [12] Höhle, M. (2020). Superspreading and the Gini Coefficient. Blog. Stockholm University, Sweden. $<$https://www.r-bloggers.com/2020/05/superspreading-and-the-gini-coefficient/$>$
  * [13] Korean, J. (1998). Distributions and moments for estimators of Gini index in an exponential distribution. _Comput. and Appl. Math._ 5(1), 213–222
  * [14] Miller, D. et.al (2020). Full genome viral sequences inform patterns of SARS-CoV-2 spread into and within Israel. _Nature Communications volume_ 11 (1)
  * [15] Pozderac, C. and Skinner, B. (2021). Superspreading of SARS-CoV-2 in the USA. _PLoS ONE_ 16(3)
  * [16] Ramasubban, T. A. (1958). The Mean Difference and the Mean Deviation of Some Discontinuous Distributions. _Biometrika_ 45, 549–556.
  * [17] Schechtman, E. and Yitzhaki, S. (2013). The Gini Methodology. _Springer Series in Statistics_ 272,481-498
  * [18] Tufekci, Z., (2020). This Overlooked Variable Is the Key to the Pandemic. _The Atlantic_ $<$https://www.theatlantic.com/health/archive/2020/09/k-overlooked-variable-driving-pandemic/616548/$>$
  * [19] White, T.J., (2007). Sharing resources: The global distribution of the Ecological Footprint. _Ecological Economics_ 64(2), 402–410.
  * [20] Yitzhaki, S. (1997). More Than a Dozen Alternative Ways of Spelling Gini. _Research on Economic Inequality._ 8, 13–30.