# Linear Regression for Power Law Distribution Fitting

Samuel<EMAIL_ADDRESS>
(December 2023)

###### Abstract

We fit the exponent of the Pareto distribution, that is equivalent or can
approximate the continuous power law distribution given a cutoff point, using
linear regression (LR). We use LR on the logged variables of the empirical
tail (one minus the empirical cumulative distribution function). We find the
distribution of the consistent LR estimator and an approximate sigmoid
relationship of the mean that underestimates the exponent. By factoring out a
sigmoid function used to approximate the mean we transform the LR estimator so
it is approximately unbiased with variance comparable to the minimum variance
unbiased transformed MLE estimator.

##  1 Introduction

Power laws in the probability distribution appear across a wide number of
disciplines [3]. We shall focus in this paper on continuous power law
distributions that approximate variables such as income, wealth and stock
returns [5]. Several definitions of a power law distribution appear in the
literature. We focus on the piecewise power law distribution and the regularly
varying distribution which we define in the next section. The power law part
of the former distribution can be exactly viewed as a Pareto distribution
whereas the power law part of the latter can only be approximated as a Pareto
distribution. We fit only the exponent of the distribution and not the cutoff
point which can be estimated using a Monte Carlo method found in Clauset et
al. [3]. There are several methods of fitting the exponent of a Pareto
distribution that are consistent (the estimator tends to the true parameter as
the sample size tends to infinity) including the method of moments, maximum
likelihood estimation (MLE), quantile methods and linear regression (LR) [11].
Many have pointed out, for example [3], that LR is inferior to MLE but to the
author’s knowledge it has not been analysed in detail why this is the case. We
focus in this paper on fitting the exponent with LR on the empirical tail (one
minus the empirical cumulative distribution funciton). We give a theoretical
result on the distribution of the LR estimator as well as approximations on
the mean and variance via simulation222Simulations are found at the author’s
github [4]. We find evidence that the LR estimator is biased with a mean that
underestimates the true exponent in a sigmoidal fashion depending on the
sample size. Factoring out a sigmoid function we find this non-linear
transformation of the LR estimator is roughly unbiased and comparable though
greater in variance to a non-linear transformation of the MLE estimator which
is unbiased of minimum variance [13]. We emphasise that the MLE estimator is
superior however we present these results on the LR estimator as they may be
of interest.

##  2 Definitions of a Power Law Distribution

As in Clauset et al. [3] we define the power law in two ways: one as a
piecewise distribution and the other asymptotically as a regularly varying
distribution. The piecewise power law distribution approximates the regularly
varying distribution after a cutoff point which we assume is known. In this
paper we will fit only the exponent of the piecewise power law distribution
assuming the cutoff. We show below how this is identical to fitting the
exponent of the Pareto distribution.

###  2.1 Piecewise Distribution

Let $X$ be a continuous random variable. We define the power law distribution
as one that has a power law tail333We note Clauset et al. [3] defines the
power law such that $\mathbb{P}(X>x)=\alpha/x^{\beta-1}$ so the density is
$f(x)=c/x^{\beta}$ for an appropriate constant $c\,.$

$\mathbb{P}(X>x)=\alpha/x^{\beta}\,,\quad x>x_{m}\,$ (1)

with $\alpha,\,\beta,\,x_{m}>0\,.$ We call $\beta$ the power law exponent and
$x_{m}$ the cutoff. Before $x_{m}$ we assume the tail of $X$ is defined by
another function so that $X$ is a piecewise distribution. By taking the
negative of the derivative of (1) we have that the density of $X$ is

$f(x)=\beta\alpha/x^{\beta+1}\,,\quad x>x_{m}\,.$ (2)

Let us define the truncated distribution:

$X_{\text{P}}=\\{X|X>x_{m}\\}\,,$ (3)

then one can integrate over the entire domain, $[x_{m},\infty)\,,$ to find
that $X_{\text{P}}$ is a Pareto distribution with tail (1) and density (2)
such that $\alpha=x_{m}^{\beta}\,.$ Now, for all $x>x_{m}\,,$

$\displaystyle\mathbb{P}(X>x)$
$\displaystyle=\mathbb{P}(X>x,X>x_{m})=\mathbb{P}(X>x_{m})\mathbb{P}(X>x|X>x_{m})$
$\displaystyle=\mathbb{P}(X>x_{m})\mathbb{P}(X_{\text{P}}>x)$
$\displaystyle=\mathbb{P}(X>x_{m})\left(\frac{x_{m}}{x}\right)^{\beta}\,.$ (4)

Thus in general $\alpha=\mathbb{P}(X>x_{m})x_{m}^{\beta}\,.$ We shall assume
the cutoff $x_{m}$ is known (or has been fitted) and we are concerned only
with fitting the exponent $\beta\,.$ We refer the reader to Clauset et al. [3]
for a Monte-Carlo method using the Kolmogorov-Smirnov statistic to fit
$x_{m}\,.$

###  2.2 Regularly Varying Distribution

We shall briefly discuss a more general asymptotic representation of the power
law distribution, the regularly varying distribution, and show how the
piecewise distribution above approximates this distribution for large enough
$x_{m}\,.$ We mention here that the method for fitting $x_{m}$ in Clauset et
al. [3] has useful statistical properties such as consistency when applied to
regularly varying distributions, see [2].

######  Definition 1 (Regularly Varying Distribution, see e.g. [7]).

A continuous random variable $X$ has a regularly varying distribution if for
all $\lambda>0$ there exists a $\beta>0$ such that the tail distribution is a
regularly varying function:

$\lim_{x\rightarrow\infty}\frac{\mathbb{P}(X>\lambda
x)}{\mathbb{P}(X>x)}=\lambda^{-\beta}\,.$

We see simply that the piecewise distribution (1) is a regularly varying
distribution.

######  Definition 2 (Slowly Varying Function, see e.g. [7, 15]).

A function $l$ is slowly varying if for all $\lambda>0$

$\lim_{x\rightarrow\infty}\frac{l(\lambda x)}{l(x)}=1\,.$

One can show (p2 of [15]) that a regularly varying distribution can be written
as

$\mathbb{P}(X>x)=l(x)x^{-\beta}$

where $l(x)$ is a slowly varying function. We have the following properties
(p18 of [15]) of slowly varying functions: for all $\varepsilon>0$

$\lim_{x\rightarrow\infty}l(x)x^{-\varepsilon}=0\,,\quad\lim_{x\rightarrow\infty}l(x)x^{\varepsilon}=\infty\,.$
(5)

We now show that a regularly varying distribution can be approximated by the
piecewise power law distribution (1).

######  Proposition 1.

If $X$ is a regularly varying distribution then for any $\varepsilon>0$ there
exists an $x_{m}>0$ such that with $\alpha=\mathbb{P}(X>x_{m})x_{m}^{\beta}$
and $x>x_{m}$

$\alpha x^{-\beta-\varepsilon}<\mathbb{P}(X>x)<\alpha
x^{-\beta+\varepsilon}\,.$

###### Proof.

Suppose for any $\varepsilon>0\,,$ that there exists an $x_{m}>0$ such that
for $x>x_{m}$ and $\alpha=\mathbb{P}(X>x_{m})x_{m}^{\beta}$ that

$\displaystyle\alpha x^{-\beta-\varepsilon}$
$\displaystyle<\mathbb{P}(X>x)<\alpha x^{-\beta+\varepsilon}\Leftrightarrow$
$\displaystyle\alpha x^{-\beta-\varepsilon}$
$\displaystyle<l(x)x^{-\beta}<\alpha x^{-\beta+\varepsilon}\Leftrightarrow$
$\displaystyle\alpha x^{-\varepsilon}$ $\displaystyle<l(x)<\alpha
x^{\varepsilon}$

which is true by the definition of the above limits of slowly varying
functions (5). ∎

Thus Proposition 1 shows that for $\varepsilon\gtrapprox 0$ there exists an
$x_{m}$ such that for $x>x_{m}$ we have $\mathbb{P}(X>x)\approx\alpha
x^{-\beta}$ with $\alpha=\mathbb{P}(X>x_{m})x_{m}^{\beta}\,.$

###  2.3 Examples

We now give an example of each of the two types of power law defined above.
First a power law distribution that is piecewise (and regularly varying) that
has an exponential then power law shape:

######  Example 1 (Piecewise power law, similar to (3.10) in [3]).

For $x_{m}\,,\,\beta>0$ define the continuous random variable $X$ by the
piecewise tail

$\mathbb{P}(X>x)=\begin{cases}e^{-\beta\frac{x}{x_{m}}}\,,&0<x\leq x_{m}\\\
e^{-\beta}\left(\dfrac{x_{m}}{x}\right)^{\beta}\,,&x>x_{m}\end{cases}$ (6)

Second a power law distribution that is not piecewise but is regularly
varying:

######  Example 2 (Lomax [9] or equivalently Pareto Type ii@ [1]).

For $\lambda\,,\,\beta>0$ define the continuous random variable $X$ by the
regularly varying tail

$\mathbb{P}(X>x)=(1+x/\lambda)^{-\beta}\,,\quad x>0\,.$ (7)

We now sample once from each of the distributions in the above examples and
plot the empirical tail, see the appendix. We fit with the unbiased minimum
variance MLE and roughly unbiased LR power law exponent estimators found in
Sections 3.1 and 3.2. For the piecewise distribution the power law part can be
viewed exactly as points from a Pareto distribution. For the Lomax
distribution a cutoff point $x_{m}$ must be chosen so that beyond $x_{m}$ we
can only approximate by a Pareto distribution. We choose the cutoff exactly
for the piecewise distribution and only approximately using the technique in
Clauset et al. [3] for the Lomax distribution. We see the fits with both
methods for the exponent of this particular Lomax distribution are not very
accurate even for a large sample size. Thus we urge caution when using these
methods for fitting the exponent of a regularly varying distribution (we find
evidence from simulation not presented here that in this case a larger initial
sample with a subsequently larger $x_{m}$ is required to get a closer estimate
consistently).

Figure 1: We generate a single sample of size $10^{4}$ from each distribution
in Example 6 and 7. For both examples we set $\beta=1.5$ and fit the power law
points (1) to estimate $\beta\,$. For Example 6 we set $x_{m}=20$ and fit the
power law, of which there are roughly $2100$ points, with both MLE (10), fit
$\hat{\beta}_{n}^{\text{MLE}}\approx 1.56$, and LR (25), fit
$\hat{\beta}_{n}^{\text{LR}}\approx 1.53\,$. For Example 7 we set $\lambda=20$
and estimate $x_{m}\approx 56$ using the method in [3] and fit the power law,
of which there are roughly $1300$ points, with both MLE (10), fit
$\hat{\beta}_{n}^{\text{MLE}}\approx 1.29$ and LR (25), fit
$\hat{\beta}_{n}^{\text{LR}}\approx 1.31\,$. Left figure is the full
distribution, right figure is where we cut off the distribution at $x_{m}$ so
that this new distribution $X_{\text{P}}$ (3) is exactly or approximated by a
Pareto distribution.

##  3 Fitting the Power Law Exponent

We consider an i.i.d. sample $\mathbf{x}=\\{x_{1},x_{2},\dots,x_{n}\\}$ from a
random variable $X$ that is from the power law tail part of the distribution
(1). We assume $x_{m}$ and $\mathbb{P}(X>x_{m})$ are known or have been fitted
and emphasise $x_{m}$ is not from the sample but is the cutoff in (1) for
which $x_{i}>x_{m}$ for all $x_{i}\in\mathbf{x}\,$. We note from (4) we only
have to fit the Pareto distribution and thus take $\mathbf{x}$ from the Pareto
distribution $X_{\text{P}}$ (3).

###  3.1 Maximum Likelihood Estimation

Suppose $X$ follows a power law distribution defined by (1). We saw that
$X_{\text{P}}$ follows a Pareto distribution with density $f(x|\beta)=\beta
x_{m}^{\beta}x^{-(\beta+1)}\,.$ Assume
$\mathbf{x}=\\{x_{1},x_{2},\dots,x_{n}\\}$ is an i.i.d. sample from
$X_{\text{P}}$ then the log-likelihood is defined

$\log L(\mathbf{x}|\beta)=\log\left(\prod_{i=1}^{n}f(x_{i}|\beta)\right)\,.$

By maximising the log-likelihood (respectively the likelihood) with respect to
$\beta$ we find that the MLE estimator for $\beta$ given $x_{m}$ is

$\hat{\beta}_{n}^{\text{MLE}_{1}}=\frac{1}{1/n\sum_{i=1}^{n}\log(x_{i}/x_{m})}\,.$
(8)

The MLE estimator (8) for the Pareto distribution has been known for a long
time with Clauset et al. [3] referencing Muniruzzaman, 1957 [10] as an early
paper containing this result. We refer the reader to [3] for a summary of
useful statistical properties of this estimator which includes consistency. We
note that (8) is exactly the Hill’s estimator for the Pareto distribution
$X_{\text{P}}$ [6]. In particular it can be shown, see Rytgaard [13], that
$\hat{\beta}_{n}^{\text{MLE}_{1}}$ has an inverse Gamma distribution with mean
and variance

$\mathbb{E}[\hat{\beta}_{n}^{\text{MLE}_{1}}]=n\beta/(n-1)\,,\quad\text{Var}(\hat{\beta}_{n}^{\text{MLE}_{1}})=n^{2}\beta^{2}/((n-1)^{2}(n-2))\,.$
(9)

This MLE estimator (8), as noted in [13], can be corrected for bias. From (9)
we see the unbiased MLE estimator is (for $n>1$)

$\tilde{\beta}_{n}^{\text{MLE}_{2}}=\frac{n-1}{n}\hat{\beta}_{n}^{\text{MLE}_{2}}=\frac{1}{\frac{1}{n-1}\sum_{i=1}^{n}\log(x_{i}/x_{m})}$
(10)

with mean and variance (for $n>2$)

$\mathbb{E}[\hat{\beta}_{n}^{\text{MLE}_{2}}]=\beta\,,\quad\text{Var}(\hat{\beta}_{n}^{\text{MLE}_{2}})=\frac{\beta^{2}}{n-2}<\text{Var}(\hat{\beta}_{n}^{\text{MLE}_{1}})\,.$
(11)

It is shown in [13] that the unbiased MLE estimator (10) has minimum variance
and is asymptotically normally distributed with mean $\beta$ and variance
$\beta^{2}/n$ and is a consistent estimator (the same for the untransformed
MLE estimator (8) [3]).

Now approximating either MLE estimator with the asymptotic normal we have that
$\approx 68$ % of values of the estimator lie within one standard deviation
i.e. lie within the interval

$(\beta-\beta/\sqrt{n},\beta+\beta/\sqrt{n})\,.$ (12)

Thus for the size of this interval to be less than $\varepsilon>0$ we want

$\frac{2\beta}{\sqrt{n}}<\varepsilon\Leftrightarrow
n>4\beta^{2}/\varepsilon^{2}\,.$

Therefore for a majority of the values of the MLE estimator to be within for
example $0.1$ of either side of $\beta$ we would want a sample of size at
least $400\beta^{2}\,.$ Now for many real-world situations $0.5<\beta<3\,,$
see e.g. Table 6.1 in [3], thus the sample is required to be quite large
before one is reasonably confident that the minimum variance MLE estimator is
close to $\beta\,.$

###  3.2 Linear Regression

We shall now apply linear regression to fit the power law exponent to the
power law model (1) using the empirical tail (one minus the empirical
distribution function, see the appendix). Let us again assume
$\mathbf{x}=\\{x_{1},x_{2},\dots,x_{n}\\}$ is an i.i.d. sample from
$X_{\text{P}}$ defined above. Then for any $x_{i}\in\mathbf{x}\,$ the
empirical tail is

$\hat{\mathbb{P}}_{n}(X_{\text{P}}\geq
x_{i})=(x_{m}/x_{i})^{\beta}+\varepsilon_{n,x_{i}}\,$ (13)

where the error $\varepsilon_{n,x_{i}}$ has the following properties (see the
appendix):

$\varepsilon_{n,x_{i}}\sim(1/n)\text{Binom}_{k\geq
1}(n,(x_{m}/x_{i})^{\beta})-(x_{m}/x_{i})^{\beta}\,,$ (14)

where $\text{Binom}_{k\geq 1}(.)$ is the truncated Binomial distribution that
restricts the support to $\\{1,2,\dots,n\\}$, for large $n$

$\varepsilon_{n,x_{i}}\approx\mathcal{N}(0,(1/n)(x_{m}/x_{i})^{\beta}(1-(x_{m}/x_{i})^{\beta}))\,$

and

$\varepsilon_{n,x_{i}}\rightarrow 0\quad\text{almost surely as
}n\rightarrow\infty\,.$ (15)

Applying log’s to (13) we have

$\displaystyle\log\hat{\mathbb{P}}_{n}(X_{\text{P}}\geq x_{i})$
$\displaystyle=\beta\log x_{m}-\beta\log
x_{i}+\log(1+(x_{m}/x_{i})^{-\beta}\varepsilon_{n,x_{i}})$
$\displaystyle=-\beta\log(x_{i}/x_{m})+\tilde{\varepsilon}_{n,x_{i}}$ (16)

Where the new error is

$\tilde{\varepsilon}_{n,x_{i}}=\log(1+(x_{m}/x_{i})^{-\beta}\varepsilon_{n,x_{i}})\,.$
(17)

We note from (14) that

$\tilde{\varepsilon}_{n,x_{i}}\sim\log B_{n,x_{i}}$ (18)

where $B_{n,x_{i}}\sim(1/n)(x_{m}/x_{i})^{-\beta}\text{Binom}_{k\geq
1}(n,(x_{m}/x_{i})^{\beta})\,.$

Applying ordinary least squares, see e.g. Section 3.1 of [14], to (16) the LR
estimator for $\beta$ can be found to be

$\hat{\beta}^{\text{OLS}_{1}}_{n}=\frac{\sum_{i=1}^{n}\log(x_{i}/x_{m})\log(1/\hat{\mathbb{P}}_{n}(X_{\text{P}}\geq
x_{i}))}{\sum_{i=1}^{n}(\log(x_{i}/x_{m}))^{2}}\,.$ (19)

We note that if $a\,,b\in\mathbb{R}_{>0}$ then by the change of base formula
for logs that

$\log_{a}x=\log_{b}x/\log_{b}a$

and so due to cancellation it does not matter which base we take for the log
in (19). As the $\tilde{\varepsilon}_{n,x_{i}}$ are a function of the
independent $x_{i}$ then they are independent. However as we will see from
simulation, the mean of (19) is biased, indicating that the errors
$\tilde{\varepsilon}_{n,x_{i}}$ do not satisfy the properties of zero mean and
constant variance to satisfy the further Gauss-Markov conditions, see e.g.
Theorem 3.2 of [14]. Due to the relatively complicated distribution of
$\tilde{\varepsilon}_{n,x_{i}}$ (18) it is however hard to prove these
properties formally. However as $n\rightarrow\infty$ we have by (15) and (17)
that

$\tilde{\varepsilon}_{n,x_{i}}\rightarrow 0\quad\text{almost surely as
}n\rightarrow\infty\,,$

thus the estimator will tend almost surely to $\beta$ as $n\rightarrow\infty$
and so the LR estimator (19) is a consistent estimator as noted by Quandt
[11].

Now let $\mathbf{x}_{\text{o}}=\\{x_{(1)},x_{(2)},\dots,x_{(n)}\\}=\mathbf{x}$
be the ordered sample of $\mathbf{x}$ thus
$x_{(1)}<x_{(2)}<\dots<x_{(n)}\,.$ We have that the empirical tail, see the
appendix, of the ordered sample is

$\hat{\mathbb{P}}_{n}(X\geq x_{(i)})=\frac{n-i+1}{n}\,,\quad i=1,2,\dots,n\,.$
(20)

As $x_{i}$ is from a Pareto distribution with parameters $x_{m}$ and $\beta$
then $k_{i}=\log(x_{i}/x_{m})$ (specifically for the natural logarithm which
we use from now on) is exponentially distributed with parameter $\beta\,$. Now
the result by Rényi, see [12], on the distribution of the order statistics for
an exponential distribution gives

$k_{(i)}=\log(x_{(i)}/x_{m})\sim\frac{1}{\beta}\sum_{j=1}^{i}\frac{Z_{j}}{n-j+1}\,,\quad
i=1,2,\dots,n\,$ (21)

where $Z_{j}\sim\text{Exp}(1)\,$ (exponentially distributed with parameter
$1\,$). Substitution of (20) and (21) into the the OLS estimator (19) leads to

$\hat{\beta}_{n}^{\text{OLS}_{1}}\sim\left(\frac{\sum_{i=1}^{n}\sum_{j=1}^{i}\frac{Z_{j}}{n-j+1}\log\left(\frac{n}{n-i+1}\right)}{\sum_{i=1}^{n}\left(\sum_{j=1}^{i}\frac{Z_{j}}{n-j+1}\right)^{2}}\right)\beta\,.$
(22)

Thus we have a representation of the distribution of the LR estimator (19). It
is a rather complicated formula which is beyond the scope of this paper to
analyse fully so we make do with analysing the LR estimator via simulation. We
note though that the term in the brackets of (22) does not depend on
$\beta\,,$ thus the distribution of the LR estimator is a random factor
multiplied by $\beta\,.$ From simulations, see Figure 2, we hypothesise two
results:

$\displaystyle\mathbb{E}[\hat{\beta}_{n}^{\text{OLS}_{1}}]$
$\displaystyle\approx\log\left(e-(\log
n)^{\gamma}/n\right)\beta\,,\quad\gamma\approx 1.6$ (23)
$\displaystyle\text{Var}(\hat{\beta}_{n}^{\text{OLS}_{1}})$
$\displaystyle\approx\mathcal{O}(1/n)$ (24)

where $e=2.718\dots$ is Euler’s number. The free parameter $\gamma$ in (23)
was approximated using non-linear least squares. Let

$r_{n}=\log\left(e-(\log n)^{\gamma}/n\right)\,,\quad\gamma\approx 1.6\,.$

This function we propose is an approximation to the mean of the random factor
in the brackets in (22). We see $r_{n}<1$ for finite $n$ and $r_{n}\rightarrow
1$ as $n\rightarrow\infty$ so that (23) underestimates $\beta$ but tends to
$\beta$ as $n\rightarrow\infty\,.$ We note that the mean approximation (23) is
sigmoidal but imperfect especially for small $n$ and further simulation would
be needed to check if it is still a good approximation above $n=1000\,,$ see
Figure 2. It would be future work to determine more about the form of (22) and
whether it has an analytic mean or variance. If an analytic mean exists then a
strictly unbiased LR estimator can be found.

Figure 2: Sample mean, left and sample variance, right of the biased LR
estimator $\hat{\beta}_{n}^{\text{OLS}_{1}}$ (19) for a Pareto distribution
$X_{\text{P}}$ with $\beta=1.5$ and known $x_{m}=1$ of 1000 samples of size
$n=10,20,\dots,1000\,.$ Black dashed line in left is the mean approximation
(23) and black dashed line in right is $\beta^{2}/n$ the asymptotic minimum
variance of the MLE estimator.

We can though still approximately correct for the bias to obtain the
approximately unbiased, consistent LR estimator

$\hat{\beta}_{n}^{\text{OLS}_{2}}=\frac{\hat{\beta}_{n}^{\text{OLS}_{1}}}{r_{n}}$
(25)

with

$\mathbb{E}[\hat{\beta}_{n}^{\text{OLS}_{2}}]\approx\beta\,,\quad\text{Var}(\hat{\beta}_{n}^{\text{OLS}_{2}})\approx\mathcal{O}(1/n)\,.$
(26)

We note though that this transformation is greater in terms of variance:

$\text{Var}(\hat{\beta}_{n}^{\text{OLS}_{2}})=\frac{\text{Var}(\hat{\beta}_{n}^{\text{OLS}_{1}})}{r_{n}^{2}}>\text{Var}(\hat{\beta}_{n}^{\text{OLS}_{1}})\,.$

However we find through simulation, see Figure 2 and Figure 3 in the next
section, that the variance of both LR estimators are not too much larger than
the minimum variance of the unbiased MLE estimator.

###  3.3 Comparison Between LR and MLE with Simulations

In Figure 3 we compare the approximately unbiased LR estimator (25) with the
minimum variance unbiased MLE estimator (10) on 1000 samples with sample sizes
of $10,20,\dots,1000.$ We show in Figure 3 that although the transformed LR
estimator is worse than the transformed MLE estimator it is comparable.

Figure 3: The approximately unbiased LR estimator
$\hat{\beta}^{\text{LR}_{2}}_{n}$ (25) and the unbiased minimum variance MLE
estimator $\hat{\beta}^{\text{MLE}_{2}}_{n}$ (10) for $\beta$ were found for a
$1000$ samples of size $n=10,20,\dots,1000$ from a Pareto distribution
$X_{\text{P}}$ with known $x_{m}=1$ and unknown power law exponent $\beta\,.$
In this case we are estimating $\beta=1.5\,.$ Top left shows the estimator
sample means against $\beta$, top right shows the sample means within one
sample standard deviation and bounded by the interval (12), bottom left shows
the sample variances where the dotted black line is the transformed MLE
estimator variance (11) of $\mathcal{O}(1/n)$ and bottom right shows
approximately the probability the LR estimator is closer to $\beta$ than the
transformed MLE estimator which fluctuates roughly around $0.42\,$.

##  4 Conclusion

We use linear regression on the empirical tail to fit the exponent of the
Pareto distribution that approximates or is identical to a continuous power
law distribution given a cutoff value. We give an analytical representation of
the distribution of the LR estimator and show evidence it’s mean
underestimates the true exponent in a sigmoidal fashion. By factoring out a
sigmoid function of the LR estimator we find comparable results to the
transformed MLE estimator. Although the minimum variance unbiased transformed
MLE estimator is superior we nonetheless present these novel results for the
LR estimator.

## References

  * Arnold, [2008] Arnold, B. C. (2008).  Pareto and generalized Pareto distributions.  In Modeling income distributions and Lorenz curves, pages 119–145. Springer.
  * Bhattacharya et al., [2020] Bhattacharya, A., Chen, B., van der Hofstad, R., and Zwart, B. (2020).  Consistency of the plfit estimator for power-law data.  arXiv preprint arXiv:2002.06870.
  * Clauset et al., [2009] Clauset, A., Shalizi, C. R., and Newman, M. E. (2009).  Power-law distributions in empirical data.  SIAM review, 51(4):661–703.
  * Forbes, [2023] Forbes, S. (2023).  GitHub simulations code.  https://github.com/saf92/power_law_fitting.
  * Gabaix, [2009] Gabaix, X. (2009).  Power laws in economics and finance.  Annu. Rev. Econ., 1(1):255–294.
  * Hill, [1975] Hill, B. M. (1975).  A simple general approach to inference about the tail of a distribution.  The annals of statistics, pages 1163–1174.
  * Jessen and Mikosch, [2006] Jessen, H. A. and Mikosch, T. (2006).  Regularly varying functions.  Publications de L’institut Mathematique, 80(94):171–192.
  * Johnson et al., [2005] Johnson, N. L., Kemp, A. W., and Kotz, S. (2005).  Univariate discrete distributions.  John Wiley & Sons.
  * Lomax, [1954] Lomax, K. S. (1954).  Business failures: Another example of the analysis of failure data.  Journal of the American statistical association, 49(268):847–852.
  * Muniruzzaman, [1957] Muniruzzaman, A. (1957).  On measures of location and dispersion and tests of hypotheses in a pare to population.  Calcutta Statistical Association Bulletin, 7(3):115–123.
  * Quandt, [1964] Quandt, R. E. (1964).  Old and mew methods of estimation and the Pareto distribution.  Metrika, 10.
  * Rényi, [1953] Rényi, A. (1953).  On the theory of order statistics.  Acta Math. Acad. Sci. Hung, 4:191–231.
  * Rytgaard, [1990] Rytgaard, M. (1990).  Estimation in the Pareto distribution.  ASTIN Bulletin: The Journal of the IAA, 20(2):201–216.
  * Seber and Lee, [2003] Seber, G. A. and Lee, A. J. (2003).  Linear regression analysis.  John Wiley & Sons.
  * Seneta, [2006] Seneta, E. (2006).  Regularly varying functions.  Springer.
  * Van der Vaart, [1998] Van der Vaart, A. W. (1998).  Asymptotic statistics.  Cambridge university press.

## Appendix: the Empirical Tail

Suppose $\mathbf{x}=\\{x_{1},x_{2},\dots,x_{n}\\}$ is an i.i.d. random sample
from any random variable $X$. The empirical cumulative distribution function
is defined, see e.g. Section 19.1 of [16]

$\hat{\mathbb{P}}_{n}(X<x)=\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{x_{i}<x}$

where $\mathbf{1}_{A}=\begin{cases}1\quad\text{ if }x\in A\\\ 0\quad\text{ if
}x\notin A\end{cases}$ is the indicator function. We define the empirical tail
as

$\hat{\mathbb{P}}_{n}(X\geq
x)=1-\hat{\mathbb{P}}_{n}(X<x)=\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{x_{i}\geq
x}\,.$ (27)

which takes values $0,1/n,2/n,\dots,1\,.$ We also introduce the error of the
empirical tail $\varepsilon_{n}$

$\hat{\mathbb{P}}_{n}(X\geq x)=\mathbb{P}(X>x)+\varepsilon_{n}\,.$

We now state some results which follow from known results on the empirical
cumulative distribution function, see e.g. Section 19.1 of [16]. We have

$\hat{\mathbb{P}}_{n}(X\geq x)\sim(1/n)\text{Binom}(n,\mathbb{P}(X\geq x))\,.$

with mean and variance

$\mathbb{E}[\hat{\mathbb{P}}_{n}(X\geq x)]=\mathbb{P}(X\geq
x)\,,\quad\text{Var}(\hat{\mathbb{P}}_{n}(X\geq x))=(1/n)\mathbb{P}(X\geq
x)\mathbb{P}(X<x)\,.$

By the law of large numbers

$\hat{\mathbb{P}}_{n}(X\geq x)\rightarrow\mathbb{P}(X\geq x)\quad\text{almost
surely as}\quad n\rightarrow\infty\,.$

By the central limit theorem

$\sqrt{n}(\hat{\mathbb{P}}_{n}(X\geq x)-\mathbb{P}(X\geq
x))\rightarrow\mathcal{N}(0,\mathbb{P}(X\geq x)\mathbb{P}(X<x))\,$

in distribution.

From the above we have the following properties for $\varepsilon_{n}\,:$

$\varepsilon_{n}\sim(1/n)\text{Binom}(n,\mathbb{P}(X\geq x))-\mathbb{P}(X\geq
x)$

with mean and variance

$\mathbb{E}[\varepsilon_{n}]=0\,,\quad\text{Var}(\varepsilon_{n})=(1/n)\mathbb{P}(X\geq
x)\mathbb{P}(X<x)\,,$ $\varepsilon_{n}\rightarrow 0\quad\text{almost surely
as}\quad n\rightarrow\infty$

and for large $n$ we can approximate the error

$\varepsilon_{n}\approx\mathcal{N}(0,(1/n)\mathbb{P}(X\geq
x)\mathbb{P}(X<x))\,.$

Let $m=\min(\mathbf{x})$ and $M=\max(\mathbf{x})$ be respectively the min and
max of the sample $\mathbf{x}\,.$ Also let
$\mathbf{x}_{\text{o}}=\\{x_{(1)},x_{(2)},\dots,x_{(n)}\\}=\mathbf{x}$ be the
ordered sample of $\mathbf{x}$ so that
$m=x_{(1)}<x_{(2)}<\dots<x_{(n)}=M\,.$ Then we have that

$\hat{\mathbb{P}}_{n}(X\geq x_{(i)})=\frac{n-i+1}{n}\in[1/n,1]\,,\quad
i=1,2,\dots,n\,.$

and in particular for $x\in[m,M]$

$\hat{\mathbb{P}}_{n}(X\geq x)\in[1/n,1]\,,\quad i=1,2,\dots,n\,.$

Thus for this case the empirical tail with $x\in[m,M]$ will not be zero and so

$\hat{\mathbb{P}}_{n}(X\geq x)\sim(1/n)\text{Binom}_{k\geq
1}(n,\mathbb{P}(X\geq x))$

where $\text{Binom}_{k\geq 1}(n,\mathbb{P}(X\geq x))$ is the truncated
Binomial distribution, see e.g. Section 3.11 of [8], excluding values that
take $0$ and thus with support $\\{1,2,\dots,n\\}\,.$


# Leptophilic Dark Matter at Linear Colliders

P. S. Bhupal Dev  Department of Physics and McDonnell Center for the Space
Sciences, Washington University, St. Louis, MO 63130, USA
E-mail<EMAIL_ADDRESS>
###### Abstract

We discuss model-independent collider constraints on the effective couplings
of leptophilic dark matter (LDM), considering its production at a future
electron-positron linear collider, with both polarized and unpolarized beam
options, in the mono-photon and mono-$Z$ channels. We show that the future
collider constraints are largely complementary to the direct and indirect
detection limits on LDM, and can potentially provide the best-ever LDM
sensitivity in the low-mass regime (below $\sim$ 10 GeV).

###### keywords:

Dark Matter, Effective Field Theory, Lepton Collider

##  1 Introduction

Many of the existing experimental constraints on dark matter (DM) crucially
rely on the DM interactions with nucleons, and therefore, can be largely
weakened if the DM predominantly interacts with the Standard Model (SM)
leptons, but not quarks at tree-level. Such leptophilic DM (LDM) could arise
naturally in many beyond the Standard Model (BSM) scenarios [1, 2, 3, 4, 5, 6,
7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22], some of which
could even explain various experimental anomalies, such as the muon anomalous
magnetic moment [23], DAMA/LIBRA annual modulation [24], anomalous cosmic ray
positron excess [25, 26, 27, 28], the galactic center gamma-ray excess [29],
and XENON1T electron excess [30]. Dedicated searches for LDM in direct
detection [31, 32, 33] and beam dump [34, 35] experiments have also been
discussed.

In this proceeding based on Ref. [36], we focus on the LDM searches at lepton
colliders, which are complementary to the direct and indirect detection
searches. We adopt an effective field theory (EFT) approach, which has been
widely used in the context of collider searches for DM following the early
works of Refs. [37, 38, 39, 40, 41, 42, 43, 44, 45]. The same interactions
responsible for DM pair-annihilation in the early universe leading to their
thermal freeze-out guarantee their direct production at colliders, as long as
kinematically allowed. This will give a characteristic mono-$X$ signature,
where the large missing transverse momentum carried away by the DM pair is
balanced by a visible sector particle $X$ (which can be either a photon, jet,
$W$, $Z$, or Higgs, depending on the model) emitted from an initial,
intermediate or final state (see Refs. [46, 47] for reviews). Specifically,
the mono-jet signature has become emblematic for LHC DM searches [48, 49, 50,
51]. However, for an LDM with loop-suppressed interactions to the SM quarks,
the hadron colliders like the LHC are not expected to provide a better limit
than the existing constraints from indirect searches, such as from AMS-02 [52,
53], at least within the EFT framework with contact interactions.

On the other hand, lepton colliders provide an ideal testing ground for the
direct production of LDM and its subsequent detection via either mono-photon
[54, 55, 56, 57, 42, 58, 59, 45, 60, 61, 62, 63] or mono-$Z$ [64, 65, 66, 67]
signatures. We go beyond the existing literature and perform a comprehensive
and comparative study of both mono-photon and mono-$Z$ signatures of LDM at
future $e^{+}e^{-}$ colliders in a model-independent, EFT approach [36]. Our
analysis is generically applicable to all future $e^{+}e^{-}$ colliders, such
as the ILC [68], CLIC [69], CEPC [70] and FCC-ee [71], but for concreteness,
we have taken the $\sqrt{s}=1$ TeV ILC as our case study for numerical
simulations. We also assume the DM to be fermionic and limit ourselves to the
dimension-6 operators, but taking into consideration all possible dimension-6
operators of scalar-pseudoscalar (S-P), vector-axialvector (V-A) and tensor-
axialtensor (T-AT) type as applicable for the most general DM-electron
coupling. Within the minimal EFT approach, the only relevant degrees of
freedom in our analysis are the DM mass and an effective cut-off scale
$\Lambda$ which determines the strength of the four-Fermi operators. This
enables us to derive model-independent ILC sensitivities on LDM in the
$(m_{\chi},\Lambda)$ plane in both mono-photon and mono-$Z$ (leptonic and
hadronic) channels, after taking into account all relevant backgrounds and
systematic uncertainties. We consider both unpolarized and polarized beam
options [68, 72], and find that with the proper choice of polarizations for
the $e^{-}$ and $e^{+}$ beams (which depends on the operator type), the DM
sensitivities could be significantly enhanced.

##  2 Effective operators

Our primary assumptions are (i) the DM particle $\chi$ couples directly only
to the SM leptons but not to the quarks (hence leptophilic), and (ii) the
energy scale of the associated new physics is large compared to the collider
energies under consideration, thus allowing us to integrate out the heavy
mediators and parametrize the DM-SM interactions using effective higher-
dimensional operators. For concreteness, we assume that the DM particles are
Dirac fermions, and therefore, the leading order DM-SM interactions are the
dimension-six four-Fermi interactions, with the most-general effective
Lagrangian given by [37]

$\displaystyle\mathcal{L}_{\rm
eff}=\frac{1}{\Lambda^{2}}\sum_{j}\left(\overline{\chi}\Gamma^{j}_{\chi}\chi\right)\left(\overline{\ell}\Gamma^{j}_{\ell}\ell\right)\,,$
(1)

where $\Lambda$ is the cut-off scale for the EFT description and the index $j$
corresponds to different Lorentz structures, as shown below. Since our main
focus is on $e^{+}e^{-}$ colliders, we will just set $\ell=e$ in Eq. (1) and
assume this to be the only leading-order coupling, but our discussion below
could be easily extended to other cases, e.g. future muon colliders [73] by
setting $\ell=\mu$.

A complete set of Lorentz-invariant operators consists of scalar (S), pseudo-
scalar (P), vector (V), axial-vector (A), tensor (T) and axial-tensor (AT)
currents. We classify them as follows:

$\displaystyle\text{S-P type}:$
$\displaystyle\Gamma_{\chi}=c^{\chi}_{S}+ic^{\chi}_{P}\gamma_{5}\,,$
$\displaystyle\Gamma_{e}=c^{e}_{S}+ic^{e}_{P}\gamma_{5}\,,$
$\displaystyle\text{V-A type}:$
$\displaystyle\Gamma_{\chi}^{\mu}=\left(c^{\chi}_{V}+c^{\chi}_{A}\gamma_{5}\right)\gamma^{\mu}\,,$
$\displaystyle\Gamma_{e\mu}=\left(c^{e}_{V}+c^{e}_{A}\gamma_{5}\right)\gamma_{\mu}\,,$
$\displaystyle\text{T-AT type}:$
$\displaystyle\Gamma_{\chi}^{\mu\nu}=\left(c^{\chi}_{T}+ic^{\chi}_{AT}\gamma_{5}\right)\sigma^{\mu\nu}\,,$
$\displaystyle\Gamma_{e\mu\nu}=\sigma_{\mu\nu}\,,$ (2)

where $\sigma^{\mu\nu}=\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$ is the spin
tensor and $c_{j}^{\chi,e}$ are dimensionless, real couplings. For simplicity,
in Eq. (1) we have used a common cut-off scale $\Lambda$ for all Lorentz
structures. Furthermore, in our subsequent numerical analysis, we will
consider one type of operator at a time, by setting the corresponding
couplings $c_{j}^{\chi,e}=1$ without loss of generality and all other
couplings equal to zero, unless otherwise specified. For instance, setting
$c_{S}^{\chi}=c_{P}^{\chi}=c_{S}^{e}=c_{P}^{e}=1$ and all other couplings
equal to zero gives us the (S+P)-type operator, which we will simply refer to
as the SP-type in the following discussion. Similarly, we will denote the
$c_{V}^{\chi}=c_{A}^{\chi}=c_{V}^{e}=c_{A}^{e}=1$ case simply as the VA-type,
and $c_{T}^{\chi}=c_{AT}^{\chi}=1$ as the TAT-type for presenting our
numerical results in the $(m_{\chi},\Lambda)$ plane. For other choices of the
couplings, our results for the sensitivity on $\Lambda$ can be easily scaled
accordingly.

We will impose a theoretical limit of $\Lambda>\sqrt{s}$ for the EFT validity.
For relatively larger DM mass, we must also have $\Lambda>2m_{\chi}$ in order
to describe DM pair annihilation by the EFT. In fact, using
$\Lambda=2m_{\chi}$ induces 100% error in the EFT prediction for $s$-channel
UV completions. Therefore, we will use $\Lambda>{\rm
max}\\{\sqrt{s},3m_{\chi}\\}$ as a conservative lower bound [74] to ensure the
validity of our EFT approach.

##  3 Mono-photon channel

For the mono-photon signal $e^{+}e^{-}\to\chi\overline{\chi}\gamma$, the
$\chi$’s will contribute to the missing transverse energy at the detector. The
dominant irreducible SM background to this process comes from neutrino pair
production with an associated ISR photon, i.e.
$e^{+}e^{-}\to\nu\overline{\nu}\gamma$. Since neutrinos are practically
indistinguishable from DMs on an event-by-event basis, the majority of this
background survives the event selection cuts. However, as we will show later,
this background is highly polarization-dependent, and therefore, can be
significantly reduced by the proper choice of polarized beams, without
affecting the signal much.

Apart from the neutrino background, any SM process with a single photon in the
final state can contribute to the total background if all other visible
particles escape detected. The SM processes containing either jets or charged
particles are relatively easy to distinguish from a DM event, so their
contribution to the total background is negligible [58]. The only exception is
the _Bhabha scattering_ process associated with an extra photon (either from
initial or final state radiation), i.e., $e^{+}e^{-}\to e^{+}e^{-}\gamma$,
which has a large cross section, is polarization-independent, and can
significantly contribute to the total background whenever the final-state
electrons and positrons go undetected, e.g. along beam pipes. In our following
analysis, we consider both neutrino and radiative Bhabha backgrounds.

###  3.1 Cross-sections

The cross-sections for the mono-photon signal
$e^{+}e^{-}\to\chi\overline{\chi}\gamma$ and the radiative neutrino background
$e^{+}e^{-}\to\nu\overline{\nu}\gamma$ at $\sqrt{s}=1$ TeV ILC are estimated
using CalcHEP [75] with proper implementation of ISR and beamsstrahlung
effects, which significantly affect the width and position of the neutrino
$Z$-resonance. For this purpose, the EFT Lagrangian (1) is implemented in
FeynRules [76] to generate the CHO library required for CalcHEP. To avoid
collinear and infrared divergences, we limit the phase space in the event
generation with the following cuts on the outgoing photon energy $E_{\gamma}$
and its polar angle $\theta_{\gamma}$:

$8\text{ GeV}<E_{\gamma}<500\text{ GeV},\;\;\;\;|\cos\theta_{\gamma}|\leq
0.995\,.$ (3)

The _radiative_ Bhabha scattering events are generated using WHIZARD [77] (to
better handle the singularities) with the same set of cuts as in Eq. (3) to
the matrix element photon (i.e., excluding the ISR and beamsstrahlung
photons). Also, some additional cuts are implemented for the Bhabha process to
take care of the soft and collinear divergences:

$\displaystyle M_{e^{\pm}_{\rm in},e^{\pm}_{\rm
out}}<2m_{e},~{}M_{e^{\pm}_{\rm out},e^{\pm}_{\rm out}}<5\text{
GeV},~{}P_{T}^{\gamma}>1\text{ GeV},~{}\Delta R_{e^{\pm},\gamma}>0.2,~{}\Delta
R_{e^{\pm},e^{\pm}}>0.4.$ (4)

Figure 1: Variation of mono-photon signal cross-section with the DM mass
(left) and the cut-off scale (right) at $\sqrt{s}=1$ TeV ILC. The solid,
dashed and dotted lines are for the SP, VA and TAT-type operators
respectively. In the left panel, the red, green and blue curves respectively
correspond to different values of the cut-off scale $\Lambda=1$ TeV, 3 TeV and
5 TeV, while in the right panel, they correspond to different values of the DM
mass $m_{\chi}=100$ GeV, 250 GeV and 450 GeV.

After generating the signal and background events, we perform a fast detector
simulation of the SiD detector of ILC [78] using Delphes3 [79] with the
configuration card validated in Ref. [80]. The variations of the unpolarized
signal cross section as a function of the DM mass and the cut-off scale are
shown in Fig. 1 left and right panels respectively for all three operator
types, namely, SP (solid), VA (dashed) and TAT (dotted)-type. We find that the
cross-section is the smallest (largest) for the SP (TAT)-type operator at any
given DM mass. In the left panel, the sudden drop in the cross-section as
$m_{\chi}$ approaches $\sqrt{s}/2$ is due to phase-space suppression.
Otherwise, for smaller DM masses, the cross-section for a given operator type
and a given cut-off scale is almost independent of the DM mass. In the right
panel, we see that for a given DM mass the cross-section drops as
$\Lambda^{-4}$, as expected.

As for the background, we find that the neutrino background cross section at
$\sqrt{s}=1$ TeV is 4.8 pb, while the radiative Bhabha background is 68.4 pb
(though it is substantially reduced after the baseline selection). On the
other hand, the DM signal cross section is found to be much smaller, as shown
in Table 3.1 for a benchmark DM mass of $m_{\chi}=100$ GeV and the cut-off
scale $\Lambda=3$ TeV.

Mono-photon background and signal cross-sections with different beam
polarizations at $\sqrt{s}=1~{}\rm{TeV}$. For the signal, we have fixed
$m_{\chi}=100~{}\rm{GeV}$ and $\Lambda=3~{}\rm{TeV}$. The numbers in bold
highlight the optimal polarization choice for a given operator type.  Process
Unpolarized Polarization Polarized cross-section (fb) type cross-section (fb)
$P(e^{-},e^{+})$ $(+,+)$ $(+,-)$ $(-,+)$ $(-,-)$ $(80,0)$ 1106 1106 8506 8506
$\nu\overline{\nu}\gamma$ $4782$ $(80,20)$ 1268 963 10160 6793 $(80,30)$ 1393
860 10993 5931 $(80,0)$ 67920 67920 68867 68867 $e^{-}e^{+}\gamma$ $68439$
$(80,20)$ 67909 68386 69285 68297 $(80,30)$ 67809 68566 69502 68181 $(80,0)$
25.5 25.5 25.5 25.5 SP-type $25.5$ $(80,20)$ 29.6 21.4 21.4 29.6 $(80,30)$
31.6 19.4 19.4 31.6 $(80,0)$ 61.7 61.7 6.9 6.9 VA-type $34.3$ $(80,20)$ 49.4
74.1 5.5 8.2 $(80,30)$ 43.2 80.3 4.8 8.9 $(80,0)$ 36.5 36.5 36.5 36.5 TAT-type
$36.5$ $(80,20)$ 42.3 30.6 30.6 42.3 $(80,30)$ 45.2 27.7 27.7 45.2

###  3.2 Effect of polarization

One important advantage of lepton colliders is that the incoming beams can be
polarized. This helps to reduce the neutrino background considerably, as shown
in Table 3.1. To utilize the full advantage of the beam polarization, we
investigate the effect of different choices of polarization on the signal and
background. At the ILC, the baseline design foresees at least 80% electron
beam polarization at the interaction point, whereas the positron beam can be
polarized up to 30% for the undulator positron source (up to 60% may be
possible with the addition of a photon collimator) [68]. For comparison, we
show our results for three different nominal absolute values of polarization:
$|P(e^{-},e^{+})|=(80,0)$, (80,20) and (80,30). In each case, we can also have
four different polarization configurations, namely, ${\rm
sign}(P(e^{-}),P(e^{+}))=(+,+)$, $(+,-)$, $(-,+)$ and $(-,-)$, where $+$ and
$-$ denote the right- and left-handed helicities respectively.

In Table 3.1, we show the effect of different schemes of polarizations and
helicity orientations on the mono-photon signal and background cross-sections.
It is clear that the radiative Bhabha background remains almost unchanged. On
the other hand, electron beam polarization is very effective in reducing the
neutrino background, as a 80% _right-handed_ electron beam can reduce the
neutrino background to $23\%$ of the unpolarized case, even without any
polarization on the positron beam. The effect is further enhanced by a _left-
handed_ positron beam. We see that for $20\%$ and $30\%$ left-handed positron
beam polarization, the neutrino background is reduced to $20\%$ and $18\%$ of
its unpolarized value, respectively.

The signals are also affected to some extent by beam polarization and the
optimal helicity configuration depends on the operator type. For SP- and TAT-
type operators we see no effect of electron-beam polarization, but a 20% (30%)
_right-handed_ positron beam can enhance the signal by 16% (24%). The VA-type
signal, on the other hand, prefers the $(+,-)$ helicity configuration – the
same choice for which the neutrino background is minimized. With the
$(+80\%,-30\%)$ configuration, the VA-type signal is enhanced by a factor of
2.3, whereas the $(+80\%,+30\%)$ configuration enhances it by a modest 26%.

Overall, although the $(+80\%,-30\%)$ configuration minimizes the background
the most, looking at the different signal to background ratio, we find that
the $(+80\%,+30\%)$ configuration is the best for the SP- and TAT-type
operators. For direct comparison between the results for different operators,
we choose to work with the $(+80\%,+30\%)$ configuration democratically for
all the operator types, unless otherwise specified.

###  3.3 Cut-based analysis

Now we analyze various kinematic distributions and perform a cut-based
analysis to optimize the signal-to-background ratio. This of course depends on
the DM mass, so in Table 3.3, we list three benchmark points (BPs) with
$m_{\chi}=100$ GeV, 250 GeV and 350 GeV respectively, and present the
corresponding selection cuts optimized for each case. Here we fix $\Lambda=3$
TeV for illustration, but in the next subsection, we will vary both $m_{\chi}$
and $\Lambda$ to obtain the $3\sigma$ sensitivity limits. As for the choice of
the DM mass values, since it was seen from Figure 1 that the signal cross-
sections are barely sensitive to the DM mass up to around $100$ GeV, our BP1
essentially captures the light DM scenario. Similarly, our BP3 is chosen
moderately close to the kinematic limit of $\sqrt{s}/2$ (going too close to
$\sqrt{s}/2$ will result in cross-section values too low too low to give
sizable event counts after all the selection cuts). The BP2 is chosen for an
intermediate mass DM in between BP1 and BP2.

Mono-photon selection cuts for different BPs across all operator types.  BP1
BP2 BP3 Definition $m_{\chi}=100\text{ GeV, }$ $m_{\chi}=250\text{ GeV, }$
$m_{\chi}=350\text{ GeV, }$ $\Lambda=3\text{ TeV}$ $\Lambda=3\text{ TeV}$
$\Lambda=3\text{ TeV}$ Baseline selection $E_{\gamma}>10\text{
GeV},\;\;|\eta_{\gamma}|<2.45,\;\;P_{T}^{\rm miss}>10\text{ GeV}$ SP-type
Cut-1 $E_{\gamma}<450\text{ GeV}$ $E_{\gamma}<340\text{ GeV}$
$E_{\gamma}<250\text{ GeV}$ Cut-2 $|\eta_{\gamma}|<1.6$ Cut-3 $P_{T}^{\rm
miss}<450\text{ GeV}$ $P_{T}^{\rm miss}<340\text{ GeV}$ $P_{T}^{\rm
miss}<240\text{ GeV}$ Cut-4 $P_{T}^{\rm frac}<1.3$ Cut-5 $1.1<\Delta
R_{\gamma,{\rm MET}}<4.5$ VA-type Cut-1 $E_{\gamma}<440\text{ GeV}$
$E_{\gamma}<350\text{ GeV}$ $E_{\gamma}<250\text{ GeV}$ Cut-2
$|\eta_{\gamma}|<1.7$ Cut-3 $P_{T}^{\rm miss}<400\text{ GeV}$ $P_{T}^{\rm
miss}<340\text{ GeV}$ $P_{T}^{\rm miss}<250\text{ GeV}$ Cut-4 $P_{T}^{\rm
frac}<1.2$ Cut-5 $1.1<\Delta R_{\gamma,{\rm MET}}<4.5$ TAT-type Cut-1
$E_{\gamma}<460\text{ GeV}$ $E_{\gamma}<360\text{ GeV}$ $E_{\gamma}<230\text{
GeV}$ Cut-2 $|\eta_{\gamma}|<1.7$ Cut-3 $P_{T}^{\rm miss}<450\text{ GeV}$
$P_{T}^{\rm miss}<350\text{ GeV}$ $P_{T}^{\rm miss}<230\text{ GeV}$ Cut-4
$P_{T}^{\rm frac}<1.2$ Cut-5 $1.1<\Delta R_{\gamma,{\rm MET}}<4.4$

We define our mono-photon signals by those events that pass through the
baseline selection criteria as defined below, in addition to the cuts given in
Eq. (3):

$E_{\gamma}>10\text{ GeV},\;\;|\eta_{\gamma}|<2.45\;\text{ and }\;P_{T}^{\rm
miss}>10\text{ GeV}\,,$ (5)

where the hardest photon in an event is considered as the signal photon. For
the radiative Bhabha background, we define the selection criteria for
electrons as $P_{T,e}>10\text{ GeV, }|\eta_{e}|<2.5$, and have kept only those
events which contain no electrons (and positrons) passing these criteria,
which means they have escaped detection. After implementing these baseline
selection cuts, we find that the signal and the neutrino background are
reduced to about 60% of their original values in Table 3.1. Similarly, the
actual Bhabha-induced background relevant for our signal is found to be only
about 13% of its original value quoted in Table 3.1 after the baseline
selection cuts, taking into account only the missed electron events. To
further enhance our signal-to-background ratio, we then examine the signal
versus background distributions of some relevant kinematic variables and
devise further cuts, which are dynamic with respect to different BPs, as
summarized in Table 3.3. See Ref. [36] for details.

Even after implementing the baseline and analysis cuts 1 through 5, the
neutrino background can only be reduced to about 40% of its original value in
Table 3.1. Similarly, the radiative Bhabha background, although substantially
reduced to about 4% of its original value in Table 3.1 after the baseline
selection and analysis cuts, still remains sizable and comparable to the
neutrino background. However, an electromagnetic calorimeter in the very
forward direction of the beamline (BeamCal) [81] can further suppress the
Bhabha background to the per mille level. To properly incorporate the effect
of BeamCal, we have used the selection efficiencies obtained from a full
detector simulation performed in Ref. [61] by modeling the complete
instrumented region in a realistic way. According to this analysis, the
selection efficiency of the Bhabha background after the BeamCal veto only is
$2.7\%$, while that of the neutrino background is between 98% and 99.6%. As
for the DM signal, we expect it to be basically unaffected (just like the
neutrino background) by the BealCal veto, as it does not contain highly
energetic charged particles in the longitudinal direction.

For the polarized case, after the baseline selection cuts, the Bhabha
background remains almost same as in the unpolarized case. The neutrino
background, on the other hand, is significantly reduced in the polarized case
to about 28% of its unpolarized value. The other cut efficiencies are also
slightly better for the neutrino background in the polarized case.

As for the signals, from Table 3.1, we see that the TAT-type operator has the
largest cross section to start with, both for the unpolarized as well as for
the $(+80\%,+30\%)$ polarized cases. Even after the baseline selection and the
specialized cuts discussed above, the TAT-type signal retains the largest
efficiency among the three types. This will be reflected in our signal
significance results below.

###  3.4 Signal significance

After implementing all the cuts mentioned above, we calculate the final signal
significance for our benchmark scenarios using the definition

${\rm Sig}=\frac{S}{\sqrt{S+B+(\epsilon B)^{2}}}\,,$ (6)

where $S$ and $B$ are the number of signal and total background events
respectively for a given integrated luminosity, and $\epsilon$ is the
background systematic uncertainty. Our results are given in Table 3.4 for the
three BPs. We show the numbers for an ideal case with zero systematics and
also for a more realistic case with 1% systematics, i.e. with $\epsilon=0.01$
(in parentheses). The results are significantly weakened in the latter case
because of the relatively large background compared to the signal.

Signal significance in the mono-photon channel for the three BPs at
$\sqrt{s}=1$ TeV. The values in parenthesis correspond to 1% background
systematic uncertainty.  Operator type Signal significance for ${\cal L}_{\rm
int}=1000\,{\rm fb}^{-1}$ Unpolarized beams Polarized beams BP-1 BP-2 BP-3
BP-1 BP-2 BP-3 SP-type $8.1\;(0.6)$ $5.8\;(0.4)$ $3.5\;(0.3)$ $18.1\;(2.4)$
$13.0\;(1.7)$ $7.8\;(1.0)$ VA-type $10.9\;(0.8)$ $8.5\;(0.6)$ $5.6\;(0.4)$
$24.9\;(3.2)$ $19.4\;(2.5)$ $12.9\;(1.7)$ TAT-type $11.8\;(0.8)$ $10.8\;(0.8)$
$8.5\;(0.6)$ $26.2\;(3.5)$ $24.1\;(3.2)$ $19.2\;(2.6)$

From Table 3.4, we see that the significance enhances as we go to lower DM
mass regions, as expected because of kinematic reasons. Operator-wise we see
that TAT and VA-type operators perform better than the SP-type. We also find
substantial (around 50%) increase in significance on application of optimal
beam polarization.

Going beyond the three BPs, we now vary the DM mass and calculate the signal
significance following the same cut-based analysis procedure outlined above.
Our results for the $3\sigma$ sensitivity contours in the $(m_{\chi},\Lambda)$
plane are shown in Figure 2 for all the operator types. The solid (dashed)
contours are for the unpolarized (optimally polarized) case, and the blue
(green) contours are assuming zero (1%) background systematics. The shaded
regions are excluded by various constraints. First of all, for $\Lambda<{\rm
max}\\{\sqrt{s}/2,3m_{\chi}\\}$, our EFT framework is not valid (cf. Sec. 2).
This is shown by the navy blue-shaded regions in Fig. 2. For $\sqrt{s}=1$ TeV
as considered here, this EFT validity limit supersedes the previous LEP limit
[42].

Figure 2: $3\sigma$ sensitivity contours in the mono-photon channel for the SP
(left), VA (middle) and TAT (right)-type operators with _unpolarized_ (solid
lines) and _polarized_ (dashed lines) $e^{+}e^{-}$ beams at $\sqrt{s}=1$ TeV
center-of-mass energy and with ${\cal L}_{\rm int}=1000$ fb-1 integrated
luminosity. The blue (green) contours are assuming zero (1%) background
systematics. The various shaded regions are excluded by direct detection
(XENON1T, PANDAX-4T), indirect detection (Fermi-LAT, AMS), astrophysics
(SN1987A) and cosmology (CMB) constraints. In the shaded region below
$\Lambda={\rm max}\\{\sqrt{s}/2,3m_{\chi}\\}$, our EFT framework is not valid.
Along the dot-dashed line, the observed DM relic density is reproduced for a
thermal DM assuming only DM-electron effective coupling.

The same effective operator given in Eq. (1) also gives rise to DM scattering
with electrons $\chi e^{-}\to\chi e^{-}$. The exact analytic expressions for
these cross sections in our EFT framework can be found in Appendix C of Ref.
[82] for all the operator types. Comparing these with the experimental upper
limits on $\sigma_{\chi e}$ from dedicated direct detection experiments [31,
32], we can derive a lower limit on the cut-off scale $\Lambda$ as a function
of the DM mass $m_{\chi}$. However, the current best limit on $\sigma_{\chi
e}$ from XENON1T is at the level of ${\cal O}(10^{-39})\>{\rm cm}^{2}$ [32],
which translates into a very weak bound on $\Lambda$ and is not relevant for
our study. Even the future ambitious proposals like DARKSPHERE can only reach
up to ${\cal O}(10^{-42})\>{\rm cm}^{2}$ [83], still 5 orders of magnitude
weaker than that needed to probe a TeV-scale $\Lambda$ value.

However, more stringent limits can be derived from DM-nucleon scattering
searches. Even for an LDM as in our case, DM-nucleon couplings are necessarily
induced at loop level from photon exchange between virtual leptons and the
quarks. In fact, as shown in Ref. [37], the loop-induced DM-nucleon scattering
almost always dominates over the DM-electron scattering. The analytic
expressions for the one and two-loop DM-nucleon scattering cross sections can
be found in Ref. [37]. We have translated the experimental upper limits from
XENON1T [84] and PANDAX-4T [85] onto the $(m_{\chi},\Lambda)$ plane, as shown
by the yellow and grey-shaded regions respectively in Fig. 2. Note that these
limits are only applicable for the vector and tensor lepton currents, i.e.
$\Gamma_{\ell}=\gamma_{\mu},\ \sigma_{\mu\nu}$ in Eq. (1). For the scalar
lepton current, $\Gamma_{\ell}=1$, the one-loop DM-nucleon coupling vanishes,
and one has to go to two loops which is suppressed by $\alpha_{\rm em}^{2}$
for the S-S type coupling and $\alpha_{\rm em}^{2}v^{2}$ (where $v\sim
10^{-3}$ is the DM velocity) for the P-S type coupling. In contrast, for
pseudo-scalar and axial-vector lepton currents, i.e.
$\Gamma_{\ell}=\gamma_{5},\ \gamma_{\mu}\gamma_{5}$, the DM-nucleon coupling
vanishes to all orders. Therefore, we have not shown the XENON1T and PANDAX-4T
limits for the SP-type operator on the top left panel of Fig. 2.

The same effective operator given in Eq. (1) also enables DM annihilation into
electrons $\chi\overline{\chi}\to e^{+}e^{-}$. The exact analytic expressions
for these cross sections in our EFT framework can be found in Appendix C of
Ref. [82] for all the operator types. Using these, we calculate the thermal-
averaged cross section times relative velocity $\langle\sigma v\rangle$ which
goes as $m_{\chi}^{2}/\Lambda^{4}$ and compare it with the existing indirect
detection upper limits on $\langle\sigma v\rangle$ in the $e^{+}e^{-}$ channel
to put a lower bound on $\Lambda$ as a function of the DM mass. This is shown
in Fig. 2 by the red and brown-shaded regions respectively for the Fermi-LAT
[86] and AMS-02 [53] constraints on $\langle\sigma v\rangle$. Similar
constraints on $\langle\sigma v\rangle$ can be derived using CMB anisotropies
[86], which is shown by the cyan-shaded region in Fig. 2, assuming an $s$-wave
annihilation (for $p$-wave annihilation, the CMB bound will be much weaker).

Along the dot-dashed line in Fig. 2, the observed relic density can be
reproduced for a DM. In principle, the region to the left and above of this
line is disfavored for a thermal DM, because in this region $\langle\sigma
v\rangle$ is smaller than the observed value of $\sim(2-5)\times 10^{-26}{\rm
cm}^{3}{\rm sec}^{-1}$ (depending on the DM mass [87]), which leads to an
overabundance of DM, since $\Omega_{\chi}h^{2}\propto 1/\langle\sigma
v\rangle$. However, this problem can be circumvented by either opening up
additional leptonic annihilation channels (like $\mu^{+}\mu^{-}$,
$\tau^{+}\tau^{-}$ and $\nu\bar{\nu}$) or even going beyond the DM paradigm
and invoking e.g., the freeze-in mechanism [88]. This will not affect our main
results, since the collider phenomenology discussed here only depends on the
DM coupling to electrons.

Also shown in Fig. 2 is the supernova constraint, which excludes the magenta-
shaded region from consideration of energy-loss and optical depth criteria
from the observation of SN1987A [82]. Here we have used an average supernova
core temperature of 30 MeV. Note that the supernova bound is only applicable
for DM mass below $\sim$ 200 MeV or so, and for a certain range of $\Lambda$
values, above which the DM particles cannot be efficiently produced in the
supernova core, and below which they will no longer free-stream.

From Fig. 2, we find that in spite of a large irreducible background, the
accessible range of the cut-off scale $\Lambda$ at $\sqrt{s}=1$ TeV ILC looks
quite promising in the mono-photon channel, especially for low mass DM, where
the collider sensitivity is almost flat, whereas the existing direct and
indirect detection constraints are much weaker. This complementarity makes the
collider searches for DM very promising. With unpolarized beams, the
3$\sigma$-reach for the SP-type operator can be up to 3.9 TeV, while for the
VA and TAT-type operators, it can be up to 4.2 TeV. With optimally polarized
beams, i.e. with $(+80\%,+30\%)$ for the SP and TAT-types and $(+80\%,-30\%)$
for the VA type, the sensitivity reaches can be extended to $4.8$ TeV (SP),
$6.5$ TeV (VA) and $5.3$ TeV (TAT), as shown in Fig. 2.

##  4 Mono-$Z$ channel

In addition to the mono-photon channel discussed in the previous section,
another useful channel for LDM search at lepton colliders is the mono-$Z$
channel, where the $Z$-boson is emitted from one of the initial states.
Depending on the subsequent decay of the $Z$-boson to either leptonic or
hadronic final states, we perform a dedicated cut-based signal and background
analysis, as discussed below.

###  4.1 Leptonic mode

For the leptonic decay of the $Z$-boson, we examine the process
$e^{+}e^{-}\to\chi\overline{\chi}Z(\to\ell^{-}\ell^{+})$. We will only
consider $\ell=e,\mu$ for simplicity and use the lepton pair as the visible
particles for tagging. The main SM background for this channel is
$e^{+}e^{-}\to\nu\overline{\nu}\;\ell^{+}\ell^{-}$, and it is polarization-
dependent.

####  4.1.1 Unpolarized and polarized cross-sections

For the signal and background simulation, we generated the UFO library for our
EFT framework using FeynRules [76] and then generated events for both signal
and background using MadGraph 5 [89] with the following basic baseline cuts:

$P_{T}(\ell)>10\text{ GeV},\;\;\;\;|\eta_{\ell}|\leq 2.5,\;\;\;\;\Delta
R_{\ell\ell}\geq 0.4\,.$ (7)

For the signal, the $Z$-bosons are decayed into the charged lepton pairs via
the MadSpin [90, 91] package which is implemented in MadGraph 5, to take care
of the spin-correlation effects of the lepton pairs. A fast detector
simulation to these events is done using Delphes 3 [79] with the same
configuration card [80] as in Sec. 3.1.

With unpolarized beams, we find that the neutrino background cross section at
$\sqrt{s}=1$ TeV is 420.5 fb, whereas the DM signal cross section is much
smaller, as shown in Table 4.1.1 for a benchmark DM mass of $m_{\chi}=100$ GeV
and the cut-off scale $\Lambda=3$ TeV. Similar to the mono-photon case, we
also examine the effect of polarization on the signal and background cross-
sections, as shown in Table 4.1.1. The neutrino background can be reduced to
28% of its original value by making the electron beam $+80\%$ polarized, and
further reduced to 21% of its original value by additionally making the
positron beam $-30\%$ polarized. The $(+80\%,-30\%)$ polarization
configuration also enhances the VA-type signal by a factor of 2.4. However,
the $(+80\%,+30\%)$ configuration is better for the SP and TAT-type signals.
For ease of comparison between different operator types, we choose to work
with the $(+80\%,+30\%)$ configuration democratically for all operator types,
as well as for the background, unless otherwise specified.

Comparison of the leptonic mono-$Z$ background and signal cross-sections for
different choices of beam polarization for $m_{\chi}=100\text{ GeV}$ and
$\Lambda=3\text{ TeV}$ at $\sqrt{s}=1\text{ TeV}$ ILC. The numbers in bold
highlight the optimal polarization choice for a given operator type.  Process
Unpolarized Polarization Polarized cross-section (fb) type cross-section (fb)
$P(e^{-},e^{+})$ $(+,+)$ $(+,-)$ $(-,+)$ $(-,-)$
$\nu\overline{\nu}\ell^{-}\ell^{+}$ $(80,0)$ $116$ $116$ $723$ $723$ $420$
$(80,20)$ $135$ $98$ $856$ $590$ $(80,30)$ $145$ $88$ $926$ $523$ SP-Type (80,
0) $0.26$ $0.26$ $0.25$ $0.25$ $0.28$ (80,20) $0.29$ $0.22$ $0.21$ $0.29$
$(80,30)$ $0.32$ $0.19$ $0.19$ $0.32$ VA-Type (80, 0) $0.15$ $0.15$ $0.02$
$0.02$ $0.08$ (80,20) $0.12$ $0.18$ $0.01$ $0.02$ $(80,30)$ $0.11$ $0.19$
$0.01$ $0.02$ TAT-Type (80, 0) $0.62$ $0.62$ $0.62$ $0.62$ $0.68$ (80,20)
$0.72$ $0.52$ $0.52$ $0.72$ $(80,30)$ $0.77$ $0.47$ $0.47$ $0.77$

####  4.1.2 Cut-based analysis

We define our signals by those events that pass through the baseline selection
criteria as defined below: $P_{T,\ell}>20\text{ GeV, }\;|\eta_{\ell}|<2.45$,
where the $Z$-boson is reconstructed by the condition that all final state
lepton-pairs are oppositely charged and of same flavor (OSSF). Other selection
criteria are dynamic with respect to different BPs. We have taken the same
three BPs as in the mono-photon case to probe different regions of the
parameter space, namely, BP1 essentially represents all light DM region, BP3
represents the region close to the kinematic limit of $\sqrt{s}/2-m_{Z}$,
whereas BP2 captures the intermediate DM mass region.

After implementing the baseline selection cuts, we find that the background is
reduced to about 40% of its original value in Table 4.1.1 for the unpolarized
(polarized) case, whereas the signals are reduced to about 60%-70% of their
original values. We then consider various kinematic distributions for the
signal and background, and devise some specialized selection cuts [36]. We
find that after applying all these cuts, we can still retain about 35%-45% of
the signal, whereas the background is reduced to below percent level of the
original values given in Table 4.1.1.

####  4.1.3 Results

After implementing all these cuts, we calculate the final signal significance
for the three BPs using Eq. (6). Our results are given in Table 4.1.3 for an
integrated luminosity of ${\cal L}_{\rm int}=1000$ fb-1. We see that as we go
higher up in the DM mass the signal significance drops. We also find that the
best-performing operator type is the TAT-type, for which more than $97\%$ of
the background events are removed after all the selection cuts. For the signal
we retain $58\%-61\%$ of the events, although for BP3 only $48\%$ remains. The
SP-type operator also gives good results, where we retain $50\%-66\%$ of the
signal across BPs and polarization choices, while removing more than $96\%$ of
the background events. Even for VA-type we retain more than $50\%$ of the
signal events and are able to cut down the background event yields to $11\%$.
We also notice the positive effect of the beam polarization by which we
achieve an enhancement of signal significance by more than $2$ times compared
to the ones with unpolarized beams. For VA-type though the significance can be
further increased for the polarized beam case by choosing the left-handed
positron beams as is evident from Table 4.1.1.

Signal significance in the mono-$Z$ leptonic channel at $\sqrt{s}=1$ TeV and
$\mathcal{L}_{\rm int}=1000\,{\rm fb}^{-1}$. The values in the parenthesis
correspond to 1% background systematic uncertainty.  Operator Type Signal
significance for ${\cal L}_{\rm int}=1000\,{\rm fb}^{-1}$ Unpolarized beams
Polarized beams BP-1 BP-2 BP-3 BP-1 BP-2 BP-3 SP-type $1.7\;(1.3)$
$0.7\;(0.5)$ $0.1\;(0.1)$ $3.8\;(3.6)$ $1.7\;(1.5)$ $0.3\;(0.3)$ VA-type
$0.2\;(0.1)$ $0.1\;(0.1)$ $0.1\;(0.1)$ $0.5\;(0.4)$ $0.4\;(0.3)$ $0.2\;(0.2)$
TAT-type $4.5\;(3.9)$ $2.4\;(1.9)$ $0.6\;(0.5)$ $9.4\;(9.1)$ $5.4\;(5.1)$
$1.3\;(1.3)$

Figure 3: $3\sigma$ sensitivity contours in the mono-$Z$ leptonic channel.
Labels are same as in Fig. 2.

Going beyond the three BPs, we now vary the DM mass and present the $3\sigma$
sensitivity reach for this channel in Fig. 3 for all the operators. The labels
and shaded regions are the same as in the mono-photon case (cf. Fig. 2). We
see that the accessible range of the cut-off scale $\Lambda$ for the
unpolarized beams can reach up to $3.2$ TeV for the TAT-type operator, whereas
for the SP and VA-type, it can reach up to 2.6 TeV and 1.6 TeV respectively.
But with the application of optimally polarized beams as discussed earlier, we
see an increase by about $25\%$ of the $3\sigma$ reach on the $\Lambda$ scale,
up to 3.2 TeV, 2.5 TeV and 4 TeV for the for SP, VA and TAT-type operators,
respectively.

###  4.2 Hadronic mode

Next we study $e^{+}e^{-}\to\chi\overline{\chi}Z(\to jj)$, where $j\equiv
u,d,c,s,b$ quarks. The relevant SM background processes for this channel are
$e^{+}e^{-}\to\nu\overline{\nu}jj$ and $e^{+}e^{-}\to jj\ell\nu$ (with one
charged lepton escaping the detector) where the jets and leptons in the final
state can come from any possible source (not necessarily from an on-shell
$Z$).

####  4.2.1 Unpolarized and polarized cross-sections

We use the same UFO library as before which is implemented using FeynRules
[76] and simulate the events for the signal and backgrounds via MadGraph 5
[89] with the following basic cuts to the parameter space:

$P_{T}(j,\ell)>10\text{ GeV},\;\;\;\;|\eta_{j}|\leq 3.0,\;\;|\eta_{\ell}|\leq
2.5,\;\;\Delta R_{jj,\ell j}\geq 0.4\,.$ (8)

For the signals, as in the leptonic case, the on-shell $Z$-bosons are decayed
into the pairs of jets using the MadSpin package [90, 91], implemented in
MadGraph 5. Both the signal and background samples are hadronized using
Pythia8.2 [92] and then the final state jets are reconstructed with with
anti-$kT$ [93] clustering algorithm with a minimum $P_{T}$ of $10$ GeV and a
cone radius ($R)$ of 0.4 using the FastJet [94] package. The fast detector
simulation to these events are done using Delphes 3 [79] with the same
configuration card [80] as discussed in Sec. 3.1.

With unpolarized beams, we find the neutrino-pair background is 798 fb,
whereas the $jj\ell\nu$ background is 1186 fb. On the other hand, the DM
signal is only at a few fb level, as shown in Table 4.2.1 for a benchmark DM
mass of $m_{\chi}=100$ GeV and the cut-off scale $\Lambda=3$ TeV. We then
examine different choices of beam-polarization on both the event samples for
this channel, as shown in Table 4.2.1. We find that both backgrounds are
polarization-dependent and fall off significantly for right-handed electron
beam and with increasing degree of polarization. We choose the polarization
configuration $P(e^{-},e^{+})=(+80\%,+30\%)$ democratically over all the
operator types.

Hadronic mono-$Z$ background and signal cross-sections for different choices
of beam polarization with $m_{\chi}=100~{}\rm{GeV}$ and $\Lambda=3~{}\rm{TeV}$
at $\sqrt{s}=1~{}\rm{TeV}$ ILC. The numbers in bold highlight the optimal
polarization choice for a given operator type.  Process Unpolarized Pol.
Polarized cross-section (fb) type cross-section (fb) $P(e^{-},e^{+})$ $(+,+)$
$(+,-)$ $(-,+)$ $(-,-)$ $(80,0)$ 178 178 1415 1415 $\nu\overline{\nu}jj$ $798$
$(80,20)$ 206 151 1689 1134 $(80,30)$ 219 136 1833 989 $(80,0)$ 302 302 2061
2061 $jj\ell\nu$ $1186$ $(80,20)$ 359 246 2446 1685 $(80,30)$ 386 216 2635
1492 $(80,0)$ 2.57 2.57 2.58 2.58 SP-Type $2.78$ $(80,20)$ 2.98 2.15 2.15 2.97
$(80,30)$ 3.17 1.95 1.95 3.17 $(80,0)$ 1.35 1.35 0.15 0.15 VA-Type $0.83$
$(80,20)$ 1.07 1.61 0.12 0.18 $(80,30)$ 0.94 1.76 0.10 0.19 $(80,0)$ 6.22 6.22
6.21 6.21 TAT-Type $6.77$ $(80,20)$ 7.23 5.23 5.24 7.21 $(80,30)$ 7.72 4.73
4.73 7.69

####  4.2.2 Cut-based analysis

After obtaining the signal and background cross-sections as reported in Table
4.2.1, we proceed with our cut-based analysis to optimize the signal
significance. We select the events that contain at least two jets with the
following transverse momentum and pseudorapidity requirements:
$P_{T,j}>20\text{ GeV, }\;|\eta_{j}|<2.45$. The hardest two jets are required
to reconstruct the $Z$-boson. Further selection cuts are applied some of which
depend on the DM mass. So, as in the leptonic channel, we have taken the same
three BPs with varying DM mass and impose dynamic cuts [36].

Signal significances of the mono-$Z$ hadronic channel at $\sqrt{s}=1~{}{\rm
TeV}$ and $\mathcal{L}_{\rm int}=1000~{}{\rm fb}^{-1}$. The values in
parenthesis correspond to 1% background systematic uncertainty.  Operator
types Signal significance for ${\cal L}_{\rm int}=1000\,{\rm fb}^{-1}$
Unpolarized beams Polarized Beam BP-1 BP-2 BP-3 BP-1 BP-2 BP-3 SP-type
$4.7\;(3.3)$ $1.5\;(1.0)$ $0.3\;(0.2)$ $10.3\;(9.1)$ $3.6\;(3.1)$ $0.8\;(0.6)$
VA-type $0.4\;(0.2)$ $0.3\;(0.1)$ $0.1\;(0.1)$ $0.2\;(0.1)$ $0.1\;(0.1)$
$0.1\;(0.04)$ TAT-type $14.2\;(10.4)$ $5.8\;(3.7)$ $1.2\;(0.7)$ $27.7\;(25.6)$
$12.9\;(11.1)$ $2.8\;(2.3)$

Figure 4: $3\sigma$ sensitivity contours in the mono-$Z$ hadronic channel.
Labels are same as in Fig. 2.

####  4.2.3 Results

The signal significances calculated using Eq. (6) are tabulated in Table
4.2.2. We see similar behavior for the different BPs as in the previously
discussed channels, i.e. enhanced signal significance with decreasing mass of
the DM. The selection cuts are most efficient for SP- and TAT-type operators.
For BP-$1$ we remove more than $98\%$ of the background events while keeping
at least $21\%$ of the signal events for the two operator types, yielding a
large signal significance especially for the TAT-type operator with polarized
beams.

Varying the DM mass, we display the $3\sigma$ sensitivity contours for all
three operators in Fig. 4. It is clear that the TAT-type operator has the best
sensitivity, which reaches up to $4.2$ TeV with unpolarized beams and $5.2$
TeV with optimally polarized beams. The SP-type operator also has a
sensitivity comparable to the mono-photon channel, and can reach up to 3.4
(4.2) TeV with unpolarized (polarized) beams. The VA-type operator has a
modest sensitivity in this channel, only up to 1.7 (2.7) TeV with unpolarized
(polarized) beams.

##  5 Conclusion

We have explored the physics potential of the future $e^{+}e^{-}$ colliders in
probing such leptophilic DM in a model-independent way. As a case study, we
have taken the $\sqrt{s}=1$ TeV ILC with an integrated luminosity of 1000 fb-1
and have analyzed the pair-production of fermionic DM using leptophilic
dimension-6 operators of all possible bilinear structures, namely, scalar-
pseudoscalar, vector-axialvector and tensor-axialtensor. We have performed a
detailed cut-based analysis for each of these operators in three different
channels based on the tagged particle, namely, mono-photon, mono-$Z$ leptonic
and hadronic.

We have taken into account one of the most important and powerful features of
lepton colliders, i.e., the possibility of beam-polarization with different
degrees of polarization and helicity orientations. We find that the ${\rm
sign}(P(e^{-}),P(e^{+}))=(+,+)$ beam configuration is optimal for the SP and
TAT-type operators, while the $(+,-)$ configuration is better for probing the
VA-type operators. The maximum value of the cut-off scale $\Lambda$ that can
be probed in each channel at $3\sigma$ is given in Table 5. We find that
without any systematics, the mono-photon channel provides the best sensitivity
across all operator types, while in presence of background systematic effects,
the mono-$Z$ hadronic channel provides better sensitivity for the SP and TAT-
type operators.

We also demonstrate the complementarity of our lepton collider study with
other existing direct and indirect detection searches for LDM (cf. Figures 2,
3 and 4). In particular, we show that lepton colliders will be able to provide
the best-ever sensitivity in the still unexplored light DM regime.

Summary of our results for the $3\sigma$ sensitivity reach of the cut-off
scale $\Lambda$ in the three different channels discussed in the text. Here we
have fixed the DM mass at 1 GeV. The numbers in parentheses are with 1%
background systematics. The numbers in bold show the highest $\Lambda$ value
that can be probed for a given operator.  Process Beam $3\sigma$ sensitivity
reach of $\Lambda$ (TeV) type configuration SP VA TAT Mono-$\gamma$
Unpolarized 3.91 (1.99) 4.19 (2.14) 4.25 (2.17) Polarized 4.84 (2.81) 6.49
(3.94) 5.28 (3.08) Mono-$Z$ Unpolarized 2.62 (2.42) 1.57 (1.36) 3.22 (3.00)
leptonic Polarized 3.24 (3.16) 2.47 (2.39) 4.02 (3.93) Mono-$Z$ Unpolarized
3.38 (2.79) 1.74 (1.39) 4.22 (3.38) hadronic Polarized 4.21 (3.87) 2.75 (2.57)
5.25 (4.71)

## Acknowledgment

The work of BD is supported in part by the US DOE Grant #DE-SC0017987.

## References

  * [1] L. M. Krauss, S. Nasri and M. Trodden, _A Model for neutrino masses and dark matter_ , _Phys. Rev. D_ 67 (2003) 085002 [hep-ph/0210389].
  * [2] E. A. Baltz and L. Bergstrom, _Detection of leptonic dark matter_ , _Phys. Rev. D_ 67 (2003) 043516 [hep-ph/0211325].
  * [3] E. Ma, _Verifiable radiative seesaw mechanism of neutrino mass and dark matter_ , _Phys. Rev. D_ 73 (2006) 077301 [hep-ph/0601225].
  * [4] T. Hambye, K. Kannike, E. Ma and M. Raidal, _Emanations of Dark Matter: Muon Anomalous Magnetic Moment, Radiative Neutrino Mass, and Novel Leptogenesis at the TeV Scale_ , _Phys. Rev. D_ 75 (2007) 095003 [hep-ph/0609228].
  * [5] R. Bernabei et al., _Investigating electron interacting dark matter_ , _Phys. Rev. D_ 77 (2008) 023506 [0712.0562].
  * [6] M. Cirelli, M. Kadastik, M. Raidal and A. Strumia, _Model-independent implications of the e+-, anti-proton cosmic ray spectra on properties of Dark Matter_ , _Nucl. Phys. B_ 813 (2009) 1 [0809.2409].
  * [7] C.-R. Chen and F. Takahashi, _Cosmic rays from Leptonic Dark Matter_ , _JCAP_ 02 (2009) 004 [0810.4110].
  * [8] X.-J. Bi, P.-H. Gu, T. Li and X. Zhang, _ATIC and PAMELA Results on Cosmic e+- Excesses and Neutrino Masses_ , _JHEP_ 04 (2009) 103 [0901.0176].
  * [9] A. Ibarra, A. Ringwald, D. Tran and C. Weniger, _Cosmic Rays from Leptophilic Dark Matter Decay via Kinetic Mixing_ , _JCAP_ 08 (2009) 017 [0903.3625].
  * [10] P. S. B. Dev, D. K. Ghosh, N. Okada and I. Saha, _Neutrino Mass and Dark Matter in light of recent AMS-02 results_ , _Phys. Rev. D_ 89 (2014) 095001 [1307.6204].
  * [11] S. Chang, R. Edezhath, J. Hutchinson and M. Luty, _Leptophilic Effective WIMPs_ , _Phys. Rev. D_ 90 (2014) 015011 [1402.7358].
  * [12] P. Agrawal, Z. Chacko and C. B. Verhaaren, _Leptophilic Dark Matter and the Anomalous Magnetic Moment of the Muon_ , _JHEP_ 08 (2014) 147 [1402.7369].
  * [13] N. F. Bell, Y. Cai, R. K. Leane and A. D. Medina, _Leptophilic dark matter with $Z^{\prime}$ interactions_, _Phys. Rev. D_ 90 (2014) 035027 [1407.3001].
  * [14] A. Freitas and S. Westhoff, _Leptophilic Dark Matter in Lepton Interactions at LEP and ILC_ , _JHEP_ 10 (2014) 116 [1408.1959].
  * [15] Q.-H. Cao, C.-R. Chen and T. Gong, _Leptophilic dark matter confronts AMS-02 cosmic-ray positron flux_ , _Chin. J. Phys._ 55 (2017) 10 [1409.7317].
  * [16] B.-Q. Lu and H.-S. Zong, _Leptophilic dark matter in Galactic Center excess_ , _Phys. Rev. D_ 93 (2016) 083504.
  * [17] G. H. Duan, L. Feng, F. Wang, L. Wu, J. M. Yang and R. Zheng, _Simplified TeV leptophilic dark matter in light of DAMPE data_ , _JHEP_ 02 (2018) 107 [1711.11012].
  * [18] E. Madge and P. Schwaller, _Leptophilic dark matter from gauged lepton number: Phenomenology and gravitational wave signatures_ , _JHEP_ 02 (2019) 048 [1809.09110].
  * [19] S. Junius, L. Lopez-Honorez and A. Mariotti, _A feeble window on leptophilic dark matter_ , _JHEP_ 07 (2019) 136 [1904.07513].
  * [20] S. Ghosh, A. Dutta Banik, E. J. Chun and D. Majumdar, _Leptophilic-portal Dark Matter in the Light of AMS-02 positron excess_ , 2003.07675.
  * [21] S. Chakraborti and R. Islam, _Implications of dark sector mixing on leptophilic scalar dark matter_ , _JHEP_ 03 (2021) 032 [2007.13719].
  * [22] S.-I. Horigome, T. Katayose, S. Matsumoto and I. Saha, _Leptophilic fermion WIMP: Role of future lepton colliders_ , _Phys. Rev. D_ 104 (2021) 055001 [2102.08645].
  * [23] Muon g-2 collaboration, _Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm_ , _Phys. Rev. Lett._ 126 (2021) 141801 [2104.03281].
  * [24] DAMA/LIBRA collaboration, _The DAMA project: Achievements, implications and perspectives_ , _Prog. Part. Nucl. Phys._ 114 (2020) 103810.
  * [25] Fermi-LAT collaboration, _Cosmic-ray electron-positron spectrum from 7 GeV to 2 TeV with the Fermi Large Area Telescope_ , _Phys. Rev. D_ 95 (2017) 082007 [1704.07195].
  * [26] DAMPE collaboration, _Direct detection of a break in the teraelectronvolt cosmic-ray spectrum of electrons and positrons_ , _Nature_ 552 (2017) 63 [1711.10981].
  * [27] CALET collaboration, _Extended Measurement of the Cosmic-Ray Electron and Positron Spectrum from 11 GeV to 4.8 TeV with the Calorimetric Electron Telescope on the International Space Station_ , _Phys. Rev. Lett._ 120 (2018) 261102 [1806.09728].
  * [28] AMS collaboration, _The Alpha Magnetic Spectrometer (AMS) on the international space station: Part II — Results from the first seven years_ , _Phys. Rept._ 894 (2021) 1.
  * [29] Fermi-LAT collaboration, _Fermi-LAT Observations of High-Energy $\gamma$-Ray Emission Toward the Galactic Center_, _Astrophys. J._ 819 (2016) 44 [1511.02938].
  * [30] XENON collaboration, _Excess electronic recoil events in XENON1T_ , _Phys. Rev. D_ 102 (2020) 072004 [2006.09721].
  * [31] XENON100 collaboration, _Exclusion of Leptophilic Dark Matter Models using XENON100 Electronic Recoil Data_ , _Science_ 349 (2015) 851 [1507.07747].
  * [32] XENON collaboration, _Light Dark Matter Search with Ionization Signals in XENON1T_ , _Phys. Rev. Lett._ 123 (2019) 251801 [1907.11485].
  * [33] LZ collaboration, _Projected sensitivities of the LUX-ZEPLIN (LZ) experiment to new physics via low-energy electron recoils_ , 2102.11740.
  * [34] C.-Y. Chen, J. Kozaczuk and Y.-M. Zhong, _Exploring leptophilic dark matter with NA64- $\mu$_, _JHEP_ 10 (2018) 154 [1807.03790].
  * [35] L. Marsicano, M. Battaglieri, A. Celentano, R. De Vita and Y.-M. Zhong, _Probing Leptophilic Dark Sectors at Electron Beam-Dump Facilities_ , _Phys. Rev. D_ 98 (2018) 115022 [1812.03829].
  * [36] S. Kundu, A. Guha, P. K. Das and P. S. B. Dev, _A model-independent analysis of leptophilic dark matter at future electron-positron colliders in the mono-photon and mono-Z channels_ , 2110.06903.
  * [37] J. Kopp, V. Niro, T. Schwetz and J. Zupan, _DAMA/LIBRA and leptonically interacting Dark Matter_ , _Phys. Rev. D_ 80 (2009) 083502 [0907.3159].
  * [38] M. Beltran, D. Hooper, E. W. Kolb, Z. A. C. Krusberg and T. M. P. Tait, _Maverick dark matter at colliders_ , _JHEP_ 09 (2010) 037 [1002.4137].
  * [39] J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd, T. M. P. Tait and H.-B. Yu, _Constraints on Light Majorana dark Matter from Colliders_ , _Phys. Lett. B_ 695 (2011) 185 [1005.1286].
  * [40] Y. Bai, P. J. Fox and R. Harnik, _The Tevatron at the Frontier of Dark Matter Direct Detection_ , _JHEP_ 12 (2010) 048 [1005.3797].
  * [41] J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd, T. M. P. Tait and H.-B. Yu, _Constraints on Dark Matter from Colliders_ , _Phys. Rev. D_ 82 (2010) 116010 [1008.1783].
  * [42] P. J. Fox, R. Harnik, J. Kopp and Y. Tsai, _LEP Shines Light on Dark Matter_ , _Phys. Rev. D_ 84 (2011) 014028 [1103.0240].
  * [43] P. J. Fox, R. Harnik, J. Kopp and Y. Tsai, _Missing Energy Signatures of Dark Matter at the LHC_ , _Phys. Rev. D_ 85 (2012) 056011 [1109.4398].
  * [44] A. Rajaraman, W. Shepherd, T. M. P. Tait and A. M. Wijangco, _LHC Bounds on Interactions of Dark Matter_ , _Phys. Rev. D_ 84 (2011) 095013 [1108.1196].
  * [45] Y. J. Chae and M. Perelstein, _Dark Matter Search at a Linear Collider: Effective Operator Approach_ , _JHEP_ 05 (2013) 138 [1211.4008].
  * [46] F. Kahlhoefer, _Review of LHC Dark Matter Searches_ , _Int. J. Mod. Phys. A_ 32 (2017) 1730006 [1702.02430].
  * [47] B. Penning, _The pursuit of dark matter at colliders—an overview_ , _J. Phys. G_ 45 (2018) 063001 [1712.01391].
  * [48] CMS collaboration, _Search for dark matter, extra dimensions, and unparticles in monojet events in proton–proton collisions at $\sqrt{s}=8$ TeV_, _Eur. Phys. J. C_ 75 (2015) 235 [1408.3583].
  * [49] ATLAS collaboration, _Search for new phenomena in final states with an energetic jet and large missing transverse momentum in pp collisions at $\sqrt{s}=$8 TeV with the ATLAS detector_, _Eur. Phys. J. C_ 75 (2015) 299 [1502.01518].
  * [50] CMS collaboration, _Search for new physics in final states with an energetic jet or a hadronically decaying $W$ or $Z$ boson and transverse momentum imbalance at $\sqrt{s}=13\text{ }\text{ }\mathrm{TeV}$_, _Phys. Rev. D_ 97 (2018) 092005 [1712.02345].
  * [51] ATLAS collaboration, _Search for new phenomena in events with an energetic jet and missing transverse momentum in $pp$ collisions at $\sqrt{s}$ =13 TeV with the ATLAS detector_, _Phys. Rev. D_ 103 (2021) 112006 [2102.10874].
  * [52] L. A. Cavasonza, H. Gast, M. Krämer, M. Pellen and S. Schael, _Constraints on leptophilic dark matter from the AMS-02 experiment_ , _Astrophys. J._ 839 (2017) 36 [1612.06634].
  * [53] I. John and T. Linden, _Cosmic-Ray Positrons Strongly Constrain Leptophilic Dark Matter_ , 2107.10261.
  * [54] DELPHI collaboration, _Photon events with missing energy in e+ e- collisions at s**(1/2) = 130-GeV to 209-GeV_ , _Eur. Phys. J. C_ 38 (2005) 395 [hep-ex/0406019].
  * [55] A. Birkedal, K. Matchev and M. Perelstein, _Dark matter at colliders: A Model independent approach_ , _Phys. Rev. D_ 70 (2004) 077701 [hep-ph/0403004].
  * [56] P. J. Fox and E. Poppitz, _Leptophilic Dark Matter_ , _Phys. Rev. D_ 79 (2009) 083528 [0811.0399].
  * [57] P. Konar, K. Kong, K. T. Matchev and M. Perelstein, _Shedding Light on the Dark Sector with Direct WIMP Production_ , _New J. Phys._ 11 (2009) 105004 [0902.2000].
  * [58] C. Bartels, M. Berggren and J. List, _Characterising WIMPs at a future $e^{+}e^{-}$ Linear Collider_, _Eur. Phys. J. C_ 72 (2012) 2213 [1206.6639].
  * [59] H. Dreiner, M. Huck, M. Krämer, D. Schmeier and J. Tattersall, _Illuminating Dark Matter at the ILC_ , _Phys. Rev. D_ 87 (2013) 075015 [1211.2254].
  * [60] Z. Liu, Y.-H. Xu and Y. Zhang, _Probing dark matter particles at CEPC_ , _JHEP_ 06 (2019) 009 [1903.12114].
  * [61] M. Habermehl, M. Berggren and J. List, _WIMP Dark Matter at the International Linear Collider_ , _Phys. Rev. D_ 101 (2020) 075053 [2001.03011].
  * [62] J. Kalinowski, W. Kotlarski, K. Mekala, P. Sopicki and A. F. Zarnecki, _Sensitivity of future $e^{+}e^{-}$ colliders to processes of dark matter production with light mediator exchange_, 2107.11194.
  * [63] B. Barman, S. Bhattacharya, S. Girmohanta and S. Jahedi, _Catch ’em all: Effective Leptophilic WIMPs at the $e^{+}\,e^{-}$ Collider_, 2109.10936.
  * [64] N. Wan, M. Song, G. Li, W.-G. Ma, R.-Y. Zhang and J.-Y. Guo, _Searching for dark matter via mono- $Z$ boson production at the ILC_, _Eur. Phys. J. C_ 74 (2014) 3219 [1403.7921].
  * [65] Z.-H. Yu, X.-J. Bi, Q.-S. Yan and P.-F. Yin, _Dark matter searches in the mono- $Z$ channel at high energy $e^{+}e^{-}$ colliders_, _Phys. Rev. D_ 90 (2014) 055010 [1404.6990].
  * [66] S. Dutta, D. Sachdeva and B. Rawat, _Signals of Leptophilic Dark Matter at the ILC_ , _Eur. Phys. J. C_ 77 (2017) 639 [1704.03994].
  * [67] B. Grzadkowski, M. Iglicki, K. Mekala and A. F. Zarnecki, _Dark-matter-spin effects at future $e^{+}e^{-}$ colliders_, _JHEP_ 08 (2020) 052 [2003.06719].
  * [68] P. Bambade et al., _The International Linear Collider: A Global Project_ , 1903.01629.
  * [69] CLIC, CLICdp collaboration, _Updated baseline for a staged Compact Linear Collider_ , 1608.07537.
  * [70] CEPC Study Group collaboration, _CEPC Conceptual Design Report: Volume 2 - Physics & Detector_, 1811.10545.
  * [71] FCC collaboration, _FCC-ee: The Lepton Collider: Future Circular Collider Conceptual Design Report Volume 2_ , _Eur. Phys. J. ST_ 228 (2019) 261.
  * [72] T. Barklow, J. Brau, K. Fujii, J. Gao, J. List, N. Walker et al., _ILC Operating Scenarios_ , 1506.07830.
  * [73] J. P. Delahaye, M. Diemoz, K. Long, B. Mansoulié, N. Pastrone, L. Rivkin et al., _Muon Colliders_ , 1901.06150.
  * [74] S. Matsumoto, S. Mukhopadhyay and Y.-L. S. Tsai, _Effective Theory of WIMP Dark Matter supplemented by Simplified Models: Singlet-like Majorana fermion case_ , _Phys. Rev. D_ 94 (2016) 065034 [1604.02230].
  * [75] A. Belyaev, N. D. Christensen and A. Pukhov, _CalcHEP 3.4 for collider physics within and beyond the Standard Model_ , _Comput. Phys. Commun._ 184 (2013) 1729 [1207.6082].
  * [76] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, _FeynRules 2.0 - A complete toolbox for tree-level phenomenology_ , _Comput. Phys. Commun._ 185 (2014) 2250 [1310.1921].
  * [77] W. Kilian, T. Ohl and J. Reuter, _WHIZARD: Simulating Multi-Particle Processes at LHC and ILC_ , _Eur. Phys. J. C_ 71 (2011) 1742 [0708.4233].
  * [78] H. Abramowicz et al., _The International Linear Collider Technical Design Report - Volume 4: Detectors_ , 1306.6329.
  * [79] DELPHES 3 collaboration, _DELPHES 3, A modular framework for fast simulation of a generic collider experiment_ , _JHEP_ 02 (2014) 057 [1307.6346].
  * [80] C. T. Potter, _DSiD: a Delphes Detector for ILC Physics Studies_ , 1602.07748.
  * [81] H. Abramowicz et al., _Forward Instrumentation for ILC Detectors_ , _JINST_ 5 (2010) P12002 [1009.2433].
  * [82] A. Guha, P. S. B. Dev and P. K. Das, _Model-independent Astrophysical Constraints on Leptophilic Dark Matter in the Framework of Tsallis Statistics_ , _JCAP_ 02 (2019) 032 [1810.00399].
  * [83] L. Hamaide and C. McCabe, _Fuelling the search for light dark matter-electron scattering_ , 2110.02985.
  * [84] XENON collaboration, _Dark Matter Search Results from a One Ton-Year Exposure of XENON1T_ , _Phys. Rev. Lett._ 121 (2018) 111302 [1805.12562].
  * [85] PandaX-4T collaboration, _Dark Matter Search Results from the PandaX-4T Commissioning Run_ , 2107.13438.
  * [86] R. K. Leane, T. R. Slatyer, J. F. Beacom and K. C. Y. Ng, _GeV-scale thermal WIMPs: Not even slightly ruled out_ , _Phys. Rev. D_ 98 (2018) 023016 [1805.10305].
  * [87] G. Steigman, B. Dasgupta and J. F. Beacom, _Precise Relic WIMP Abundance and its Impact on Searches for Dark Matter Annihilation_ , _Phys. Rev. D_ 86 (2012) 023506 [1204.3622].
  * [88] L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ 03 (2010) 080 [0911.1120].
  * [89] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ , _JHEP_ 07 (2014) 079 [1405.0301].
  * [90] S. Frixione, E. Laenen, P. Motylinski and B. R. Webber, _Angular correlations of lepton pairs from vector boson and top quark decays in Monte Carlo simulations_ , _JHEP_ 04 (2007) 081 [hep-ph/0702198].
  * [91] P. Artoisenet, R. Frederix, O. Mattelaer and R. Rietkerk, _Automatic spin-entangled decays of heavy resonances in Monte Carlo simulations_ , _JHEP_ 03 (2013) 015 [1212.3460].
  * [92] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159 [1410.3012].
  * [93] M. Cacciari, G. P. Salam and G. Soyez, _The anti- $k_{t}$ jet clustering algorithm_, _JHEP_ 04 (2008) 063 [0802.1189].
  * [94] M. Cacciari, G. P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J. C_ 72 (2012) 1896 [1111.6097].