# Operator-difference schemes on non-uniform grids for second-order
evolutionary equations

P.N. Vabishchevich11footnotemark: 122footnotemark: 2<EMAIL_ADDRESS>Nuclear
Safety Institute, Russian Academy of Sciences, 52, B. Tulskaya, 115191 Moscow,
Russia  North-Eastern Federal University, 58, Belinskogo st, Yakutsk, 677000,
Russia

###### Abstract

The approximate solution of the Cauchy problem for second-order evolution
equations is performed, first of all, using three-level time approximations.
Such approximations are easily constructed and relatively uncomplicated to
investigate when using uniform time grids. When solving applied problems
numerically, we should focus on approximations with variable time steps. When
using multilevel schemes on non-uniform grids, we should maintain accuracy by
choosing appropriate approximations and ensuring the approximate solution’s
stability. In this paper, we construct unconditionally stable first- and
second-order accuracy schemes on a non-uniform time grid for the approximate
solution of the Cauchy problem for a second-order evolutionary equation. We
use a special transformation of the original second-order differential-
operator equation to a system of first-order equations. For the system of
first-order equations, we apply standard two-level time approximations. We
obtained stability estimates for the initial data and the right-hand side in
finite-dimensional Hilbert space. Eliminating auxiliary variables leads to
three-level schemes for the initial second-order evolution equation. Numerical
experiments were performed for the test problem for a one-dimensional in space
bi-parabolic equation. The accuracy and stability properties of the
constructed schemes are demonstrated on non-uniform grids with randomly
varying grid steps.

###### keywords:

Second-order evolutionary equation , Cauchy problem , System of first-order
evolutionary equations , Two-level schemes , Stability of the approximate
solution , Bi-parabolic equation

###### MSC:

35R20 , 65J08 , 65M06 , 65M12

††journal: arXiv

##  1 Introduction

Two-level schemes (Runge-Kutta methods) [1, 2] are widely used for the
approximate solution of the Cauchy problem for systems of ordinary
differential equations. The second class of methods is based on multilevel
time approximations (multistep methods) when transitioning to a new time level
involves a known solution on two or more preceding time levels. Two- and
three-level approximations are most common when numerically solving
nonstationary boundary value problems for partial differential equations. When
considering the Cauchy problem for second-order evolutionary equations, it is
natural to focus on three-level approximations.

Multilevel schemes on uniform time grids are the easiest to construct and
investigate. The computational practice focuses on using non-uniform grids for
accounting for local features of the solution by choosing the time step. When
using multilevel approximations on non-uniform grids, there are two main
problems. The first one is related to the choice of approximations in time to
ensure the required accuracy of the approximated solution. The second problem
is generated by ensuring the conditions for the stability of the solution
under the selected approximations.

The results of the theory of stability (correctness) of operator-difference
schemes [3] are widely used for the construction and investigation of time-
difference approximations in the numerical solution of applied problems. In
this case, the difference scheme is considered in finite-dimensional Hilbert
spaces; the scheme is written in a unified (canonical) form for which
stability conditions exist. The necessary and sufficient conditions for
stability of two- and three-level schemes on uniform time grids have been
obtained (see, for example, [4]).

For two-level schemes, the transition to non-uniform grids is editorial. In
the case of three-level schemes, the situation is much more complicated, both
in the choice of approximations and in the study of stability. Let us note
some available results in this direction in the particular case when the time
step change is rare. When considering the stability of three-level difference
schemes, we must consider the problem separately when the time step increases
and when it decreases. We have found (see, e.g., [5, 6, 7]) that when the time
step changes, the a priori stability estimates become worse, with the constant
in the a priori estimate increasing when the time step changes in proportion
to the time step ratio.

We can proceed to a system of first-order evolution equations for high-order
evolution equations. We can use the usual two-level approximations on non-
uniform grids for the resulting vector problem. The specificity of the problem
is considered by choosing auxiliary unknowns so that the two-level schemes for
the system of first-order equations have the necessary properties of accuracy
and stability. This paper implements such an approach for the Cauchy problem
for a general second-order differential-operator equation with self-adjoint
constant positively determined operators in a finite-dimensional Hilbert
space. Unconditionally stable first- and second-order approximation schemes
for the introduced system of first-order evolution equations have been
constructed. Algebraically equivalent unconditionally stable three-level
schemes of the second order are obtained when the auxiliary variable is
eliminated.

We give a brief content of the paper. Section 2 sets up a Cauchy problem for a
second-order evolutionary equation. A typical result on the stability of a
three-level difference scheme on a uniform time grid is established. The
original second-order equation is written down in Section 3 as a system of
first-order equations when a new unknown quantity is introduced. An a priori
estimate for the solution of the resulting vector Cauchy problem is proved.
The critical part of the paper is Section 3. The unconditionally stable two-
level schemes for the system of equations are constructed. With the
elimination of auxiliary variables, a three-level scheme on a non-uniform time
grid was obtained for the original second-order evolution equation. Section 5
presents the results of numerical experiments on a uniform and non-uniform
grids. The Cauchy problem for the one-dimensional in space bi-parabolic
equation is considered a test problem. We summarize the results of our study
in Section 6.

##  2 Problem formulation

We consider the Cauchy problem for a second-order evolutionary equation in the
real finite-dimensional Hilbert space $H$. The function $u(t)$ satisfies the
equation

$C\frac{d^{2}u}{dt^{2}}+B\frac{du}{dt}+Au=f(t),\quad t>0,$ (2.1)

and the initial conditions

$u(0)=u_{0},\quad\frac{du}{dt}(0)=u^{\prime}_{0}.$ (2.2)

The linear constants (independent of $t$) operators $A,B,C$ are self-adjoint
and positively definite:

$A=A^{*}\geq\nu_{A}I,\quad\nu_{A}>0,\quad
B=B^{*}\geq\nu_{B}I,\quad\nu_{B}>0,\quad C=C^{*}\geq\nu_{C}I,\quad\nu_{C}>0,$
(2.3)

where $I$ is a unit operator in $H$. We will use the usual notations
$(\cdot,\cdot)$ and $\|\cdot\|$ for the scalar product and norm in $H$. For a
self-adjoint positively defined operator $D$, we define a Hilbert space
$H_{D}$ with scalar product and norm $(u,v)_{D}=(Du,v),\
\|u\|_{D}=(u,v)_{D}^{1/2}.$

We get such problems, in particular, when using the finite element method or
the finite volume method for approximation by space in the numerical study of
nonstationary boundary value problems for partial differential equations. We
can associate equation (2.1) at $C=0$ with a parabolic equation and a
hyperbolic equation at $B=0$.

Let us give a simple a priori estimate for the solution of the problem
(2.1)—(2.3); this estimate will guide us in the study of time approximations.
By multiplying equation (2.1) scalarly in $H$ by $du/dt$, we obtain the
equality

$\frac{1}{2}\frac{d}{dt}\Big{(}\Big{\|}\frac{du}{dt}\Big{\|}_{C}^{2}+\|u\\\
|A^{2}\Big{)}+\Big{\|}\frac{du}{dt}\Big{\|}_{B}^{2}=\Big{(}f,\frac{du}{dt}\Big{)}.$

For the right-hand side, we use the estimate

$\Big{(}f,\frac{du}{dt}\Big{)}\leq\Big{\|}\frac{du}{dt}\Big{\|}_{B}^{2}+\frac{1}{4}\|f\\\
|2_{B^{-1}}.$

Given this inequality, Gronwall’s lemma gives us an a priori estimate

$\Big{\|}\frac{du}{dt}(t)\Big{\|}_{C}^{2}+\|u(t)\|_{A}^{2}\leq\|v_{0}\|_{C}^{2}+\|u_{0}\|_{A}^{2}+\frac{1}{2}\int_{0}^{t}\|f(\theta)\|^{2}_{B^{-1}}d\theta.$
(2.4)

In the approximate solution of the Cauchy problem (2.1), (2.2), it is natural
to focus on unconditionally stable three-level difference schemes. It is
easiest to work with uniform time grids when the time step $\tau=\mathop{\rm
const}\nolimits$ and let $y^{n}=y(t_{n}),\ t_{n}=n\tau$, $n=0,1,\ldots$. A
three-level difference scheme with weight $\sigma=\mathrm{const}>0$ has the
form

$C\frac{y^{n+1}-2y^{n}+y^{n-1}}{\tau^{2}}+B\frac{y^{n+1}-y^{n-1}}{2\tau}+Ay^{n+\sigma}=f^{n}\quad
n=1,2,\ldots,$ (2.5) $y^{0}=u^{0},\quad y^{1}=\accentset{\sim}{u}^{1},$ (2.6)

when using the notations

$y^{n+\sigma}=\sigma y^{n+1}+(1-2\sigma)y^{n}+\sigma y^{n-1}.$

Given the initial conditions (2.2), we put

$u^{0}=u_{0},\quad\accentset{\sim}{u}^{1}=u_{0}+\tau
u^{\prime}_{0}\quad(\accentset{\sim}{u}^{1}\approx u^{1}(\tau)\approx
u_{0}+\tau u^{\prime}_{0}).$

The three-level scheme (2.5), (2.6) has the second order of approximation by
$\tau$. The theory of stability of operator-difference schemes gives the
following conditions for unconditional stability of this scheme: $\sigma\geq
0.25$.

######  Theorem 1.

The three-level difference scheme (2.3), (2.5), (2.6) is unconditionally
stable for $\sigma\geq 0.25$. Under these constraints, for the approximate
solution of (2.1), (2.2), the a priori estimate is

$\begin{split}\Big{\|}\frac{y^{n+1}-y^{n}}{\tau}\Big{\|}_{D}^{2}&+\Big{\|}\frac{y^{n+1}+y^{n}}{2}\Big{\|}^{2}_{A}\\\
&\leq\Big{\|}\frac{\accentset{\sim}{u}^{1}-u^{0}}{\tau}\Big{\|}_{D}^{2}+\Big{\|}\frac{\accentset{\sim}{u}^{1}+u^{0}}{2}\Big{\|}^{2}_{A}+\frac{1}{2}\sum_{k=1}^{n}\tau\|f^{k}\|^{2}_{B^{-1}},\quad
n=1,2,\ldots,\end{split}$ (2.7)

where

$D=C+\Big{(}\sigma-\frac{1}{4}\Big{)}\tau^{2}A.$

######  Proof 1.

Given the equality

$\sigma y^{n+1}+(1-2\sigma)y^{n}+\sigma
y^{n-1}=\frac{1}{4}(y^{n+1}+2y^{n}+y^{n-1})+\Big{(}\sigma-\frac{1}{4}\Big{)}(y^{n+1}-2y^{n}+y^{n-1}),$

we write equation (2.5) as

$D\frac{y^{n+1}-2y^{n}+y^{n-1}}{\tau^{2}}+B\frac{y^{n+1}-y^{n-1}}{2\tau}+\frac{1}{4}A(y^{n+1}+2y^{n}+y^{n-1})=f^{n}.$
(2.8)

Let’s denote by

$s^{n}=y^{n}-y^{n-1},\quad r^{n}=\frac{1}{2}(y^{n}+y^{n-1}),$

then

$\begin{split}y^{n+1}-y^{n-1}&=s^{n+1}+s^{n},\\\
y^{n+1}-2y^{n}+y^{n-1}&=s^{n+1}-s^{n},\\\
y^{n+1}+2y^{n}+y^{n-1}&=2(r^{n+1}+r^{n}).\end{split}$

Equation (2.8) takes the form

$D\frac{s^{n+1}-s^{n}}{\tau^{2}}+B\frac{s^{n+1}+s^{n}}{2\tau}+A\frac{r^{n+1}+r^{n}}{2}=f^{n}.$
(2.9)

Let’s multiply equation (2.9) by $s^{n+1}+s^{n}=2(r^{n+1}-r^{n})$; this give

$\frac{1}{\tau^{2}}\big{(}\|s^{n+1}\|_{D}^{2}-\|s^{n}\|_{D}^{2}\big{)}+\|r^{n+1}\|^{2}_{A}-\|r^{n}\|^{2}_{A}+\frac{1}{2\tau}\|s^{n+1}+s^{n}\|_{B}^{2}=(f^{n},s^{n+1}+s^{n}).$

Using the introduced notations and the inequality

$(f^{n},s^{n+1}+s^{n})\leq\frac{1}{2\tau}\|s^{n+1}+s^{n}\|_{B}^{2}+\frac{\tau}{2}\|f^{n}\|^{2}_{B^{-1}},$

we get an estimate of the solution on one level in time

$\Big{\|}\frac{y^{n+1}-y^{n}}{\tau}\Big{\|}_{D}^{2}+\Big{\|}\frac{y^{n+1}+y^{n}}{2}\Big{\|}^{2}_{A}\leq\Big{\|}\frac{y^{n}-y^{n-1}}{\tau}\Big{\|}_{D}^{2}+\Big{\|}\frac{y^{n}+y^{n-1}}{2}\Big{\|}^{2}_{A}+\frac{1}{2}\tau\|f^{n}\|^{2}_{B^{-1}}.$
(2.10)

We obtain the provable estimate (2.7) by applying the difference analog of
Gronwall’s lemma.

######  Remark 1.

We can consider the a priori estimate (2.7) as a direct discrete analogue of
the estimate (2.4) for solving the Cauchy problem (2.1)–(2.2).

In computational practice, in simulating nonstationary processes, we should
focus on using non-uniform time grids. For multilevel schemes on non-uniform
grids, there are problems with maintaining the accuracy of the approximations.
The stability conditions of schemes on non-uniform grids also deteriorate. In
particular, when considering three-level operator-difference schemes, we have
stability estimates that depend on the time step ratio and how often the time
step change occurs. Our goal is to construct nonuniform time approximations
for nonstationary problems (2.1)–(2.3), which would give unconditionally
stable schemes of the second order of accuracy.

##  3 System of first-order equations

Constructing approximations on non-uniform grids for first-order evolution
equations is solved using various variants of two-level schemes (Runge-Kutta
methods). On the other hand, we can always proceed from evolution equations of
the second and higher order to a system of evolution equations of the first
order. In this case, the focus should be on choosing new variables. We
consider such possibilities of constructing approximations on non-uniform time
grids for the approximate solution of the problem (2.1)–(2.3).

The simplest version of introducing a new unknown quantity $v=du/dt$ for
equation (2.1), (2.3) is not of great interest. We put

$v=C\frac{du}{dt}+Bu$

and write equation (2.1) as a system of equations

$C\frac{du}{dt}+Bu-v=0,$ (3.1) $A^{-1}\frac{dv}{dt}+u=A^{-1}f,\quad t>0.$
(3.2)

For the system (3.1), (3.2), a Cauchy problem is posed when the initial
conditions are

$u(0)=u_{0},\quad v(0)=v_{0},$ (3.3)

where given (2.2) we have $v_{0}=Cu^{\prime}_{0}+Bu_{0}$.

It is convenient for us to write the system of equations (3.1), (3.2) as one
equation for vector quantities. Define vectors $\bm{u}=\\{u,v\\}$ and
$\bm{f}=\\{0,A^{-1}f\\}$ and from (3.1)–(3.2) we obtain the Cauchy problem

$\bm{B}\frac{d\bm{u}}{dt}+\bm{A}\bm{u}=\bm{f},\quad t>0,$ (3.4)
$\bm{u}(0)=\bm{u}_{0},$ (3.5)

where $\bm{u}_{0}=\\{u_{0},v_{0}\\}$. For the operator matrices $\bm{B}$ and
$\bm{A}$, we obtain the expressions

$\bm{B}=\mathrm{diag}\,\big{(}C,A^{-1}\big{)},\quad\bm{A}=\left(\begin{array}[]{cc}B&-I\\\
I&0\\\ \end{array}\right).$

The problem (3.4), (3.5) we will consider on the direct sum of spaces
$\bm{H}=H\oplus H$. For $\bm{u}_{1},\bm{u}_{2}\in\bm{H}$ the scalar product
and norm are defined by the expressions

$(\bm{u}_{1},\bm{u}_{2})=(u_{1},v_{1})+(u_{2},v_{2}),\quad\|\bm{u}\|=(\bm{u},\bm{u})^{1/2}.$

Given (2.3), we obtain

$\bm{B}=\bm{B}^{*}>0,\quad\bm{A}\geq 0.$ (3.6)

Such good properties of the $\bm{B}$ and $\bm{A}$ operators are provided by
choosing the second equation of the system as (3.2).

To obtain an a priori estimate for the solution of (3.4), (3.5), we multiply
the equation (3.4) by $\bm{u}$ scalarly in $\bm{H}$. Considering (3.6), this
gives

$\|\bm{u}\|_{\bm{B}}\frac{d}{dt}\|\bm{u}\|_{\bm{B}}\leq(\bm{f},\bm{u}).$

Taking into account of

$(\bm{f},\bm{u})\leq\|\bm{\varphi}\|_{\bm{B}^{-1}}\|\bm{u}\|_{\bm{B}},$

we get

$\|\bm{u}(t)\|_{\bm{B}}\leq\|\bm{u}_{0}\|_{\bm{B}}+\int_{0}^{t}\|\bm{f}(s)\|_{\bm{B}^{-1}}ds.$
(3.7)

In our case

$\|\bm{u}(t)\|_{\bm{B}}=\Big{(}\|u(t)\|_{C}^{2}+\|v(t)\|_{A^{-1}}^{2}\Big{)}^{1/2},\quad\|\bm{v}_{0}\|_{\bm{B}}=\Big{(}\|u_{0}\|_{C}^{2}+\|v_{0}\|_{A^{-1}}^{2}\Big{)}^{1/2},\quad\|\bm{f}(t)\|_{\bm{B}^{-1}}=\|f(t)\|_{A^{-1}},$

so for the problem (3.1)–(3.3) we have an estimate

$\Big{(}\|u(t)\|_{C}^{2}+\|v(t)\|_{A^{-1}}^{2}\Big{)}^{1/2}\leq\Big{(}\|u_{0}\|_{C}^{2}+\|v_{0}\|_{A^{-1}}^{2}\Big{)}^{1/2}+\int_{0}^{t}\|f(s)\|_{A^{-1}}ds.$
(3.8)

######  Remark 2.

For the primary problem (2.1)–(2.3), the estimate (3.8) is established by
scalarly multiplying equation (2.1) by

$A^{-1}v=A^{-1}\Big{(}C\frac{du}{dt}+Bu\Big{)}.$

######  Remark 3.

The positive definiteness condition of the operator $B$ (see (2.3)) can be
relaxed. The estimates (3.7), (3.8) take place even when $\nu_{B}=0$. In this
case, we must replace the estimate of (2.4) for the solution of (2.1), (2.2)
by (3.8).

##  4 Operator-difference schemes on non-uniform grids

Two-level schemes are most easily constructed to approximate the problem
(3.4), (3.5). We consider approximations on a non-uniform time grid when

$t_{n+1}=t_{n}+\tau_{n+1},\quad n=0,1,\ldots,\quad t_{0}=0.$

For the approximate solution, we will use a two-level scheme with weight
$\sigma=\mathop{\rm const}\nolimits$

$\bm{B}\frac{\bm{y}^{n+1}-\bm{y}^{n}}{\tau_{n+1}}+\bm{A}\bm{y}^{n+\sigma}=\bm{f}^{n+\sigma}\quad
n=0,1,\ldots,$ (4.1) $\bm{y}^{0}=\bm{u}_{0},$ (4.2)

when using the notations

$\bm{y}^{n}=(y^{n},w^{n}),\quad\bm{y}^{n+\sigma}=\sigma\bm{y}^{n+1}+(1-\sigma)\bm{y}^{n}.$

For a sufficiently smooth solution of $\bm{u}(t)$ of the problem (3.4), (3.5),
the difference scheme (4.1), (4.2) has first order approximation by $\tau$ for
$\sigma\neq 0.5$. For $\sigma=0.5$ we have a Crank-Nicolson scheme; this
scheme approximates (3.4), (3.5) with second order on $\tau$. Necessary and
sufficient stability conditions for such schemes are well known in the theory
of stability (correctness) of operator-difference schemes. The following
statement guarantees sufficient conditions for unconditional stability.

######  Theorem 2.

At $\sigma\geq 0.5$ the two-level operator-difference scheme (3.4), (3.5) is
unconditionally stable. In this case, there is an a priori estimate for the
approximate solution

$\|\bm{y}^{n+1}\|_{\bm{B}}\leq\|\bm{u}_{0}\|_{\bm{B}}+\sum_{k=0}^{n}\tau_{n+1}\|\bm{f}^{k+\sigma}\|_{\bm{B}^{-1}},\quad
n=0,1,\ldots.$ (4.3)

######  Proof 2.

We will multiply equation (4.1) scalarly by $\tau\bm{y}^{n+\sigma}$. Given the
non-negativity of the operator $\bm{A}$, we have the inequality

$(\bm{B}(\bm{y}^{n+1}-\bm{y}^{n}),\bm{y}^{n+\sigma})\leq\tau\|\bm{f}^{n+\sigma}\|_{\bm{B}^{-1}}\|\bm{y}^{n+\sigma}\|_{\bm{B}}.$
(4.4)

At $\sigma\geq 0.5$ we use (see lemma 1 in [8]) for the left-hand side of the
inequality (4.4) the estimate

$\big{(}\bm{B}(\bm{y}^{n+1}-\bm{y}^{n}),y^{n+\sigma}\big{)}\geq\big{(}\|\bm{y}^{n+1}\|_{\bm{B}}-\|\bm{y}^{n}\|_{\bm{B}}\big{)}\|\bm{y}^{n+\sigma}\|_{\bm{B}}.$

By doing so, we get an estimate for the solution on the new time level

$\|\bm{y}^{n+1}\|_{\bm{B}}\leq\|\bm{y}^{n}\|_{\bm{B}}+\tau\|\bm{f}^{n+\sigma}\|_{\bm{B}^{-1}}.$

The difference analog of Gronwall’s lemma gives us an a priori estimate (4.3),
which is the grid analog of the estimate (3.7).

Applying the vector scheme (3.4), (3.5) for the approximate solution of the
problem (3.1)–(3.3) leads us to the operator-difference scheme

$C\frac{y^{n+1}-y^{n}}{\tau_{n+1}}+B(\sigma y^{n+1}+(1-\sigma)y^{n})-(\sigma
w^{n+1}+(1-\sigma)w^{n})=0,$ (4.5) $\frac{w^{n+1}-w^{n}}{\tau_{n+1}}+A(\sigma
y^{n+1}+(1-\sigma)y^{n})=f^{n+\sigma},$ (4.6) $y^{0}=u_{0},\quad w^{0}=v_{0}.$
(4.7)

Under these conditions, the stability estimate (4.3) takes the form

$\Big{(}\|y^{n+1}\|_{C}^{2}+\|w^{n+1}\|_{A^{-1}}^{2}\Big{)}^{1/2}\leq\Big{(}\|u_{0}\||_{C}^{2}+\|v_{0}\|_{A^{-1}}^{2}\Big{)}^{1/2}+\sum_{k=0}^{n}\tau_{k+1}\|f^{k+\sigma}\|_{A^{-1}}.$
(4.8)

The estimate (4.8) is the discrete analog of the a priori estimate (3.8) for
the problem (3.1)–(3.3).

The computational realization of the transition to a new level in time for the
system of equations (3.1), (3.2) is provided as follows. From equation (4.6),
we get the representation

$w^{n+1}=-\sigma\tau_{n+1}Ay^{n+1}+\chi_{1}^{n},$ (4.9)

at a known

$\chi_{1}^{n}=w^{n}-(1-\sigma)\tau_{n+1}Ay^{n}+\tau_{n+1}f^{n+\sigma}.$

Similarly, from equation (4.5) we have

$(C+\sigma\tau_{n+1}B)y^{n+1}=\sigma\tau_{n+1}w^{n+1}+\chi_{2}^{n},$ (4.10)

where

$\chi_{2}^{n}=(C-(1-\sigma)\tau_{n+1}B)y^{n}+(1-\sigma)\tau_{n+1}w^{n}.$

Substituting (4.9) into (4.10) gives the equation to determine $y^{n+1}$

$(C+\sigma\tau_{n+1}B+\sigma^{2}\tau^{2}_{n+1}A)y^{n+1}=\chi^{n},$ (4.11)

at the right-hand side

$\chi^{n}=\sigma\tau_{n+1}\chi_{1}^{n}+\chi_{2}^{n}.$

After solving the equation (4.11), the expression (4.9) is used to determine
$w^{n+1}$.

When applying the standard three-level scheme on a uniform grid (2.5), (2.6)
while determining $y^{n+1}$ we solve the problem

$\Big{(}C+\frac{1}{2}\tau
B+\sigma\tau^{2}A\Big{)}y^{n+1}=\widetilde{\chi}^{n},$

with the right-hand side

$\widetilde{\chi}^{n}=C(2y^{n}-y^{n-1})+\frac{1}{2}\tau
By^{n-1}-\tau^{2}A((1-2\sigma)y^{n}+\sigma y^{n-1})+\tau^{2}f^{n}.$

A comparison with the problem (4.11) at $\tau_{n+1}=\tau$ shows that these
problems are close (left-hand sides are the same) if we put $\sigma=0.25$ in
equation (2.5) and set $\sigma=0.5$ in the system (4.5), (4.6).

The two-level scheme (4.5), (4.6) for two variables $y$ and $w$ can be written
as a three-level scheme only for $y$. To do this, we write the equation (4.5)
on two adjacent time levels and write their difference in the form

$\begin{split}C\Big{(}\frac{y^{n+1}-y^{n}}{\tau_{n+1}}-\frac{y^{n}-y^{n-1}}{\tau_{n}}\Big{)}&+B\big{(}\sigma(y^{n+1}-y^{n})+(1-\sigma)(y^{n}-y^{n-1})\big{)}\\\
&=\sigma(w^{n+1}-w^{n})+(1-\sigma)(w^{n}-w^{n-1}).\end{split}$ (4.12)

For the right-hand side, the equation (4.6) is also used at two levels in
time; this gives

$\begin{split}\sigma(w^{n+1}-w^{n})&+(1-\sigma)(w^{n}-w^{n-1})\\\
&=-\sigma\tau_{n+1}A(\sigma
y^{n+1}+(1-\sigma)y^{n})-(1-\sigma)\tau_{n}A(\sigma
y^{n}+(1-\sigma)y^{n-1})\\\
&+\sigma\tau_{n+1}f^{n+\sigma}+(1-\sigma)\tau_{n}f^{n-1+\sigma}.\end{split}$
(4.13)

From (4.12), (4.13), we get a three-level scheme for the variable $y$.

In the most interesting case $\sigma=0.5$, when the scheme (4.5), (4.6) has
the second order of accuracy, we obtain

$\begin{split}C\frac{2}{\tau_{n}+\tau_{n+1}}\Big{(}\frac{y^{n+1}-y^{n}}{\tau_{n+1}}-\frac{y^{n}-y^{n-1}}{\tau_{n}}\Big{)}&+B\Big{(}\overline{\sigma}_{n+1}\frac{y^{n+1}-y^{n}}{\tau_{n+1}}+(1-\overline{\sigma}_{n+1})\frac{y^{n}-y^{n-1}}{\tau_{n}}\Big{)}\\\
&+A\Big{(}\overline{\sigma}_{n+1}\frac{y^{n+1}+y^{n}}{2}+(1-\overline{\sigma}_{n+1})\frac{y^{n}+y^{n-1}}{2}\Big{)}\\\
&=\overline{\sigma}_{n+1}f^{n+1/2}+(1-\overline{\sigma})_{n+1}f^{n-1/2}.\end{split}$
(4.14)

The variable weight $\overline{\sigma}_{n+1}$ in (4.14) is determined by the
expression

$\overline{\sigma}_{n+1}=\frac{\tau_{n+1}}{\tau_{n}+\tau_{n+1}}.$

Given that for a sufficiently smooth function $\varphi(t)$

$\overline{\sigma}_{n+1}\varphi(t^{n+1/2})+(1-\overline{\sigma})_{n+1}\varphi(t^{n-1/2})=\varphi(\overline{t}^{n})+\mathcal{O}\big{(}(\tau_{n+1}+\tau_{n})^{2}\big{)},\quad\overline{t}^{n}=\frac{1}{2}(t^{n+1}+t^{n-1}),$

we relate the scheme (4.14) to an approximation of equation (2.1) on an non-
uniform grid at $t=\overline{t}^{n}$.

##  5 Numerical experiments

Among the generalized heat conduction models [9], we will distinguish the bi-
parabolic equation [10]. In this case, the homogeneous second-order
evolutionary equation has the form

$\Big{(}\frac{d}{dt}+D\Big{)}u+\alpha\Big{(}\frac{d}{dt}+D\Big{)}^{2}u=0,\quad
t>0,$ (5.1)

where

$D=D^{*}\geq\nu_{D}I,\quad\nu_{D}>0,\quad\alpha=\mathop{\rm
const}\nolimits>0.$

We take the initial conditions for equation (5.1) as

$u(0)=u_{0},\quad\frac{du}{dt}(0)=0.$ (5.2)

Equation (5.1) is written as a second-order evolutionary equation (2.1) at

$C=\alpha I,\quad B=I+2\alpha D,\quad A=D+\alpha D^{2}.$

The solution of the problem (5.1), (5.2) is

$u(t)=\Big{(}I+\alpha\Big{(}1-\exp\Big{(}-\frac{t}{\alpha}\Big{)}D\Big{)}\exp(-Dt)u_{0}.$
(5.3)

It follows from the representation (5.3) that the relaxation of the initial
state by the exponential law is provided not only by the positively defined
operator $D$ but also by the square of the parabolic operator (parameter
$\alpha$).

The principal goal of our study is to construct time approximations. We can
trace the key points when considering one-dimensional space problems. It is
optional to move to more interesting multidimensional problems. We use
sufficiently detailed computational grids over space when we can neglect
computational errors of approximations over space.

Let us assume that in dimensionless variables $x\in[0,1]$ and introduce a
uniform grid with step $h$. Let us denote by $\omega$ the set of internal grid
nodes in space. The test problem is considered under homogeneous boundary
conditions of the first kind. The difference operator $-D$ is associated with
the second derivative operator on $x$. For $y(x)=0,\ x\notin\omega$, we define
the operator $D$ by the relation

$Dy=-\frac{1}{h^{2}}\big{(}y(x+h)-2y(x)+y(x-h)\big{)},\quad x\in\omega.$

To be able to compare the calculation results with data from other generalized
thermal conductivity models (see [11]), we restrict ourselves to the initial
conditions

$u_{0}(x)=\left\\{\begin{array}[]{rr}x,&0<x\leq 0.5,\\\ \ 0,&0.5<x<1,\\\
\end{array}\right.\quad x\in\omega.$

We performed calculations on a grid over space with $h=2\cdot 10^{-3}$ for
$0<t\leq T,\ T=0.1$. We used the scheme (4.5), (4.6) with $\sigma=0.5$. For
the second initial condition (4.7) taking into account of (5.2), we have
$v_{0}=u_{0}+2\alpha Du_{0}$.

The new level problem (4.11) for finding $y^{n+1}$ under our conditions is

$Ry^{n+1}=\chi^{n},$

where

$R=\alpha I+\frac{\tau_{n+1}}{2}(I+2\alpha
D)+\frac{\tau^{2}_{n+1}}{4}(D+\alpha D^{2}).$

For the operator $R$, we use the factorization

$R=(pD+I)\big{(}\alpha pD+(p+\alpha)I\big{)},\quad p=\frac{\tau_{n+1}}{2}.$

To determine $y^{n+1}$, two common discrete problems with operators $D+cI$ at
corresponding values of $c=\mathop{\rm const}\nolimits>0$ are solved.

Figure 1: Exact solution of the problem at different points in time when
$\alpha=0$. Figure 2: Exact solution of the problem at different points in
time when $\alpha=0.01$. Figure 3: Exact solution of the problem at different
points in time when $\alpha=0.1$.

The exact solution of the problem (5.1), (5.2) on the used by us grid over
space at different values of the parameter $\alpha$ is shown in Figures 1–3.
The solution of the standard parabolic problem ($\alpha=0$, Fig.1) is
characterized by a monotonic decrease over time of the positive maximum of the
solution, which is a consequence of the maximum principle. The most striking
manifestation of the hyperbolicity of the problem at positive $\alpha$ is (see
Fig.2 and Fig.3) the appearance of a region of negative solutions.

The error of the approximate solution of the problem (5.1), (5.2) was
estimated in the grid space norm $L_{2}(\omega)$ for particular time moments
$t=t^{n}$ as follows:

$\varepsilon(t^{n})=\Big{(}\sum_{x\in\omega}(y^{n}(x)-u(x,t^{n}))^{2}h\Big{)}^{1/2}.$

The accuracy of the approximate solution using a uniform time grid
($\tau_{n}=\tau=T/N$) is illustrated by Fig.4. We observed the convergence to
be approximately second order when a factor of two and four increases the
number of steps $N$.

Figure 4: The error of the approximate solution of the problem (5.1), (5.2) at
$\alpha=0.01$ on different time grids.

It is convenient to illustrate the capabilities of the constructed schemes on
non-uniform time grids with the results obtained by randomly varying the grid
step. The rule calculated the time step

$\tau_{n}=\tau\big{(}1+q(\xi-0.5)\big{)},\quad
n=1,2,\ldots,N,\quad\tau=\frac{T}{N},$

where $\xi$ is a random variable uniformly distributed on the interval
$[0,1]$. In the calculations below, the parameter $q=0.5$; thus, the grid
steps can vary by $5/3$. The accuracy of the approximate solution using a non-
uniform random time grid is shown in Fig.5. The non-uniform grids used are
shown in Figure 6. Comparison of Figures 4 and 5 shows that the use of a
substantially non-uniform time grid has little effect on the accuracy of the
approximate solution of the considered Cauchy problem for the second-order
evolution equation.

Figure 5: The error of the approximate solution at $\alpha=0.01$ on different
time non-uniform grids. Figure 6: Time grid steps.

##  6 Conclusions

  1. 1.

A Cauchy problem is posed for a second-order operator evolution equation in a
finite-dimensional Hilbert space with self-adjoint constant operators. We have
well-known results about unconditional stability and convergence of schemes
with weights on uniform grids for such problems. The problem of constructing
and investigating three-level operator-difference schemes using non-uniform
grids in terms of time is considered.

  2. 2.

A second-order evolution equation is written down as a system of two evolution
equations with an appropriate choice for new variables. Two-level schemes with
weights on non-uniform time grids are used for approximate solutions to vector
problems. Conditions for unconditional stability of the solution of the Cauchy
problem for the system of equations on the initial data and the right-hand
side in the corresponding Hilbert spaces have been established. An equivalent
three-level scheme for the original second-order evolution equation has been
obtained.

  3. 3.

Numerical results of the approximate solution of the Cauchy problem for the
evolution equation on non-uniform time grids are presented. The problem for
the bi-parabolic equation, which generalizes classical models of heat
conduction, is considered a test one. The computational data on a grid with
randomly changing steps demonstrate the performance of the proposed time
approximations.

## Acknowledgements

This work has been supported by the grants the Russian Science Foundation (RSF
23-41-00037 and RSF 23-71-30013)

## References

  * [1] J. C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, Chichester, 2008.
  * [2] E. Hairer, G. Wanner, Solving Ordinary Differential Equations. II: Stiff and Differential-Algebraic Problems, Springer, Berlin, 1996.
  * [3] A. A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, New York, 2001\.
  * [4] A. A. Samarskii, P. P. Matus, P. N. Vabishchevich, Difference Schemes with Operator Factors, Kluwer Academic, Dordrecht, 2002.
  * [5] A. A. Samarskii, P. N. Vabishchevich, E. L. Makarevich, P. P. Matus, Stability of three-layer difference schemes on time-nonuniform grids, Doklady Mathematics 63 (1) (2001) 106–108.
  * [6] P. Matus, E. Zyuzina, Three-level difference schemes on non-uniform in time grids, Computational Methods in Applied Mathematics 1 (3) (2001) 265–284.
  * [7] C.-H. Cho, Stability for the finite difference schemes of the linear wave equation with nonuniform time meshes, Numerical Methods for Partial Differential Equations 29 (3) (2013) 1031–1042.
  * [8] P. N. Vabishchevich, Flux-splitting schemes for parabolic equations with mixed derivatives, Computational Mathematics and Mathematical Physics 53 (8) (2013) 1139–1152.
  * [9] D. D. Joseph, L. Preziosi, Heat waves, Reviews of Modern Physics 61 (1) (1989) 41\.
  * [10] V. I. Fushchich, A. S. Galitsyn, A. S. Polubinskii, A new mathematical model of heat conduction processes, Ukrainian Mathematical Journal 42 (1990) 210–216.
  * [11] P. N. Vabishchevich, Numerical solution of the heat conduction problem with memory, Computers & Mathematics with Applications 118 (2022) 230–236.


.eps.gz .eps.gzpdf.pdfzcat #1 | epstopdf –filter –outfile=

###### Abstract

The direct pair-production of the superpartner of the $\tau$-lepton, the
$\widetilde{\tau}$, is one of the most interesting channels to search for SUSY
in. First of all, the $\widetilde{\tau}$ is likely to be the lightest of the
scalar leptons. Secondly the signature of $\widetilde{\tau}$ pair production
is one of the experimentally most difficult ones, thereby constituting the
“worst” possible scenario for SUSY searches. The current model-independent
$\widetilde{\tau}$ limits comes from analyses performed at LEP but they suffer
from the limited energy of this facility. Limits obtained at the LHC do extend
to higher masses, but they are only valid under strong assumptions. The
International Linear Collider, the ILC, is a future electron-positron
collider, foreseen to operate initially at a centre-of-mass energy of 250 GeV,
then to be upgraded to 500 GeV, and possibly to 1 TeV at a later stage. ILC
will be a powerful facility for SUSY searches. The capability of the ILC for
determining exclusion/discovery limits for the $\widetilde{\tau}$ in a model-
independent way is shown in this paper, together with an overview of the
current state-of-the-art. Results of the last studies of $\widetilde{\tau}$
pair-production at the ILC are presented, showing the improvements with
respect to previous results. A detailed study of the “worst” scenario for
$\widetilde{\tau}$ exclusion/discovery, taking into account the effect of the
$\widetilde{\tau}$ mixing on $\widetilde{\tau}$ production cross-section and
detection efficiency, is presented. The study also includes an analysis of the
effect of the overlay particles in the $\widetilde{\tau}$ searches. The
conclusion is that both the exclusion and discovery reaches for this “worst”
case would extend to only a few GeV below the kinematic limit at the ILC. Also
scenarios with the $\widetilde{\tau}$ and the Lightest SUSY Particle (the LSP)
quite close in mass can be discovered or excluded at most $\widetilde{\tau}$
masses. The studies were done using detailed detector simulation of the ILD
concept at the ILC. For signal, the fast detector simulation SGV was used,
while the full Geant4 based DDSim was used for the standard model backgrounds.

Submitted to the Proceedings of the US Community Study

on the Future of Particle Physics (Snowmass 2021)

##  1 Introduction

Supersymmetry (SUSY) [1, 2, 3, 4, 5] is one of the most promising candidates
for new physics. It could explain or at least give some hint at solutions to
current problems of the Standard Model (SM), such as the gauge hierarchy
problem, the nature of Dark Matter (DM) or the possible theory-experiment
discrepancy of the muon magnetic moment. SUSY is a symmetry of spacetime
relating fermions and bosons. For every SM particle it introduces a
superpartner with the same quantum numbers except for the spin. The spin
differs by half a unit from the value of its SM partner. A new parity,
R-parity, is commonly introduced in SUSY, which has a crucial impact in SUSY
phenomenology. R-parity takes an even value for SM particles and odd value for
the SUSY ones. Multiplicative R-parity conservation 111The introduction and
conservation of this symmetry is inspired by flavour physics constraints since
the most general SUSY Lagrangian induces flavour-changing neutral interactions
that are avoided imposing R-parity conservation., assumed in most of the SUSY
models, implies that the SUSY particles are always created in pairs and that
the lightest SUSY particle (the LSP) is stable and, when cosmological
constraints are taken into account, also neutral. An important point in this
kind of studies is the fundamental SUSY principle stating that the couplings
of particles and sparticles are related by symmetry. This allows to know the
cross sections for SUSY pair production solely from the knowledge of initial
centre-of-mass energy of the collider and the masses of the involved SUSY
particles.

##  2 SUSY searches

Considerable efforts have been and are being devoted to the search of SUSY at
different facilities. Searches at hadron colliders, such as the LHC, are
mainly sensitive to the production of coloured sparticles, i.e. the gluino and
the squarks. They are most probably the heaviest ones. The search of the
lighter colour-neutral SUSY states, such as sleptons, charginos or
neutralinos, at hadron colliders is challenged by the much smaller cross
sections, and high backgrounds. Mass limits have been obtained at the LHC, but
they are only valid if many constraints on the model parameters are fulfilled.
Lepton colliders, like LEP, have a higher sensitivity to the production of
colour-neutral SUSY states, but the searches up to now were limited by the
beam energy. Limits computed at these facilities are however valid for any
value of the model parameters not shown in the exclusion plots. The future
International Linear Collider (ILC) [6], an electron-positron collider
operating at centre-of-mass energies of $250-500$ GeV and with upgrade
capability to $1$ TeV, is seen as an ideal environment for SUSY studies. SUSY
searches at the ILC would profit from the high electron and positron beam
polarisations, $\pm$80$\%$$\mp$30$\%$ 222We introduce the following notation
for beam polarisations,
$\mathcal{P}\equiv(\mathcal{P}_{e^{-}},\mathcal{P}_{e^{+}})$, and define the
pure beam polarisations as $\mathcal{P}_{LR}\equiv(-1,+1)$ and
$\mathcal{P}_{RL}\equiv(+1,-1)$. The nominal beam polarisations for the ILC
are defined as $\mathcal{P}_{-+}\equiv(-0.8,+0.3)$ and
$\mathcal{P}_{+-}\equiv(+0.8,-0.3)$. in this study, a well defined initial
state (in 4-momentum and spin configuration), a clean and reconstructable
final state, near absence of pile-up, hermetic detectors (almost 4$\pi$
coverage) and trigger-less operation, which is a huge advantage for precision
measurements and unexpected signatures.

##  3 Motivation for $\widetilde{\tau}$ searches

For evaluating the power of SUSY searches at future facilities, it is
beneficial to focus on the lightest particle in the SUSY spectrum that could
be accessible. Since the cosmological constraints requires a neutral and
colourless LSP, the next-to-lightest SUSY particle, the NLSP, would be the
first one to be detected. The NLSP can only decay to the LSP and the SM
partner of the NLSP (or “virtually” via its SM partner, if the LSP-NLSP mass-
difference is smaller than the mass of the SM partner). This already makes the
NLSP production special: heavier states might well decay in cascades, and thus
have signatures that depend strongly on the model. Furthermore, there is only
a finite set of sparticles that could be the NLSP, so a systematic search for
each possible case is feasible. This also means that one can a priori estimate
which will be the most difficult case, namely the NLSP that combines small
production cross-section with a difficult experimental signature. The
$\widetilde{\tau}$ satisfies both these conditions. Therefore, studies of
$\widetilde{\tau}$ production might be seen as the way to determine the
guaranteed discovery or exclusion reach for SUSY: any other NLSP would be
easier to find.

The $\widetilde{\tau}$ is the super-partner of the ${\tau}$. Like for any
other fermion or sfermion, there are two weak hyper-charge states,
$\widetilde{\tau}_{R}$ and $\widetilde{\tau}_{L}$. For the fermions the chiral
symmetry assures that both weak hyper-charge states are degenerate in mass.
However this symmetry does not apply to sfermions, since they are scalars,
rather than fermions. Hence there is no reason to expect that
$\widetilde{\tau}_{R}$ and $\widetilde{\tau}_{L}$ would have the same mass.
Furthermore, mixing between the weak hyper-charge states yields the physical
states. The strength of the couplings involved in the mixing of states depend
on the fermion mass and hence only the third generation of the sleptons,
$\widetilde{\tau}$, will mix 333This is also the case for the squarks, where
the third generation, the stop the and the sbottom, are expected to mix.. As a
consequence of the mixing the lightest $\widetilde{\tau}$,
$\widetilde{\tau}_{1}$, would most likely be the lightest slepton, due to the
seesaw mechanism: the mass of the lightest physical (mixed) state would be
smaller than the mass of any un-mixed weak hyper-charge state. The cross
section of the $\widetilde{\tau}$ also differs from the one of the ${\tau}$,
not only due to phase-space limitations - the $\widetilde{\tau}$ being more
massive than the ${\tau}$ \- but also due to the mixing. In $e^{+}e^{-}$
colliders, assuming R-parity conservation, the $\widetilde{\tau}$ will be
pair-produced, with contribution of the s-channel only, via $Z^{0}$/$\gamma$
exchange. The strength of the $Z^{0}$/$\gamma$ $\widetilde{\tau}$
$\widetilde{\tau}$ coupling depends on the $\widetilde{\tau}$ mixing, reaching
its minimum value when the coupling $\widetilde{\tau}_{1}$
$\widetilde{\tau}_{1}$ $Z^{0}$ vanishes, which will lead to the worst possible
scenario in $\widetilde{\tau}$ and, in general, in slepton searches. Another
property making the search of $\widetilde{\tau}$ the worst case, is the fact
that its SM partner is unstable. It decays before it can be detected, and, as
a further complication, some of its decay products are undetectable neutrinos.
This on one hand makes the identification more difficult than the direct decay
to electrons or muons, and on the other hand, since the decay products are
only partially detectable, that kinematic signatures get blurred. The search
of a light $\widetilde{\tau}$ is also theoretically motivated: SUSY models
with a light $\widetilde{\tau}$ can accommodate the observed relic density, by
enhancing the $\widetilde{\tau}$-neutralino co-annihilation process [7].

##  4 Limits at LEP, LHC and previous ILC studies

The most model-independent limit on the $\widetilde{\tau}$ mass comes from the
LEP experiments [8]. They set a minimum value that ranges from 87 to 93 GeV
depending on the mass difference between the $\widetilde{\tau}$ and the
neutralino, not smaller than 7 GeV. These limits, shown in figure 1(a), are
valid for any mixing and any value of the model-parameters, other than the two
masses explicitly shown in the plot.

(a) (b)

Figure 1: (a): 95$\%$ CL exclusion limits for $\widetilde{\tau}$ pair
production obtained combining data collected at the four LEP experiments with
energies ranging from 183 GeV to 208 GeV. From [8]. (b): 95$\%$ CL exclusion
and discovery potential for $\widetilde{\tau}$ pair production at the HL-LHC,
assuming $\widetilde{\tau}_{L}^{+}\widetilde{\tau}_{L}^{-}$ \+
$\widetilde{\tau}_{R}^{+}\widetilde{\tau}_{R}^{-}$ production,
$\widetilde{\tau}_{L}^{+}\widetilde{\tau}_{L}^{-}$ production or
$\widetilde{\tau}_{R}^{+}\widetilde{\tau}_{R}^{-}$ production. From [9].

An analysis by the DELPHI experiment, targeted at low mass differences,
excludes a $\widetilde{\tau}$ with mass below 26.3 GeV, for any mixing, and
any mass difference larger than the $\tau$ mass [10].

At the LHC, ATLAS and CMS have determined limits on the $\widetilde{\tau}$
mass, analysing data from Run 1 and Run 2 [11, 12]. These limits, however, are
only valid under certain assumptions. Both experiments assume
$\widetilde{\tau}_{R}$ and $\widetilde{\tau}_{L}$ to be mass-degenerate. This
is a very unlikely scenario, the running of the $\widetilde{\tau}_{R}$ and
$\widetilde{\tau}_{L}$ masses from the GUT scale to the weak scale follows
renormalisation group equations with $\beta$-functions that are inevitably
different for the two weak hyper-charge states. They also assume that there is
no mixing between the weak hyper-charge eigenstates, which is again very
improbable. The mixing will yield to cross section of the lightest physical
state smaller than that of any unmixed state. Putting together
$\widetilde{\tau}_{R}$ and $\widetilde{\tau}_{L}$ by adding the cross
sections, ATLAS excludes $\widetilde{\tau}$ masses between approximately 120
and 390 GeV for a nearly mass-less neutralino444100$\%$ decay to
$\widetilde{\tau}$ and neutralino is assumed, as it is in the analysis
presented in this paper. Under the same conditions, CMS extends the lower
limit to 90 GeV closing the gap with the LEP limit. Analysis of a pure
$\widetilde{\tau}_{L}$ pair production set limits between 150 and 310 GeV from
ATLAS data and up to 125 GeV from CMS; both limits again assume a nearly mass-
less neutralino.

The future HL-LHC should provide an improvement on the $\widetilde{\tau}$
limits provided by ATLAS and CMS, not only because of an increase of the
luminosity but also because of an expected gain in sensitivity to direct
$\widetilde{\tau}$ production due to the use of different analysis methods.
Simulation studies have already been performed in both experiments [9, 13].
Upper limits for $\widetilde{\tau}$ masses are indeed increased by about 300
GeV, but they suffer from the same constraints as the previous studies. ATLAS
adds limits for pure $\widetilde{\tau}_{R}$ pair production, that could be
considered the closest case to the physical lightest $\widetilde{\tau}$ since
it is likely to be the lightest of the two weak hyper-charge states and is the
one with the lower cross section. These limits, presented in figure 1(b), show
that no discovery potential is expected in this case, only exclusion
potential. They do not have exclusion potential for $\widetilde{\tau}$ co-
annihilation scenarios, a highly motivated scenario if SUSY is to provide a
viable DM candidate: Such a scenario requires that the $\widetilde{\tau}$-LSP
mass difference is small, $\lesssim$ 10 GeV.

$\widetilde{\tau}$ searches at the ILC have been also performed in previous
studies [14]. They assume an integrated luminosity of 500 fb-1 at
$\sqrt{s}=500$ GeV and only used the $\mathcal{P}_{+-}$ data sample. In the
current study, both $\mathcal{P}_{+-}$ and $\mathcal{P}_{-+}$ data are used,
and the luminosity is increased to the one corresponding the the foreseen
running scenario, 1.6 ab-1 at each of the beam-polarisation combinations. The
previous study was aimed at scanning the entire $M_{LSP}-M_{NLSP}$ plane, and
doing so for several different NLSP candidates. Quite generic cuts were
therefore used, and were optimised “on-the-fly” at different points. More
specifically, the limits presented in that study do not have a dedicated
analysis for low mass differences between the $\widetilde{\tau}$ and the LSP,
$\Delta M$, and are only valid down to $\Delta M$ 3-4 GeV. The exclusion limit
goes up to 240 GeV with a discovery potential up to 230 GeV for large mass
differences.

A recent study of slepton production at ILC and, in general, future e+e-
colliders can be found in [15]. They have demonstrated that these colliders
would be able to probe most of the kinematically accessible parameter space
with only a few days of data, with the capability of discovering/excluding new
physics that would evade detection at the LHC.

##  5 Conditions and tools

The study was done assuming an integrated luminosity of 1.6 ab-1 at
$\sqrt{s}=500$ GeV for each of the beam polarisations $\mathcal{P}_{+-}$ and
$\mathcal{P}_{-+}$, according to the H-20 running scenario for the ILC500
[16]555The H20 scenario is defined as: $\sqrt{s}$=500 GeV, total integrated
luminosity 4 ab-1 with 1.6 ab-1 for $\mathcal{P}_{+-}$ and $\mathcal{P}_{-+}$
, 0.4 ab-1 for $\mathcal{P}_{--}$ and $\mathcal{P}_{++}$.. The study assumes
R-parity conservation and a 100$\%$ decay of the $\widetilde{\tau}$ to
${\tau}$ and the lightest neutralino ($\chi^{0}_{1}$), the LSP in this case.
In order to select the worst scenario, the $\widetilde{\tau}$ mixing angle was
set to 53 degrees. This mixing angle corresponds to the lowest cross section
in the case of un-polarised beams, due to the suppression of the s-channel
with $Z$ exchange in the $\widetilde{\tau}$ pair production. In section 7.3,
we will show that this mixing angle in fact also corresponds to the worst case
at ILC operating according to the H20 scenario. The generated background event
samples were those of the standard “IDR” production [17]. They were generated
with Whizard 1.95 [18], and contain all standard model processes with up to
six fermions in the final state. Beam-spectra and the amount of photons in the
beams were simulated with GuineaPig [19]. Detector simulation and
reconstruction were done on the Grid using DDSim [20] and Marlin [21]. The
Grid production was done by the ILD production team using the Dirac [22]
system. The SGV fast detector simulation [23], adapted to the ILD concept [17]
at ILC, was used for detector simulation and event reconstruction for signal
events. These events were generated with Whizard v2.8.5, interfaced to Tauola
[24]. Tauola simulates the ${\tau}$-decays taking into account the ${\tau}$
polarisation in the products. The same beam-spectrum as was used for the fully
simulated background samples was also used for the signal sample. Both the
signal and background samples were analysed using the tools included in SGV,
and the relevant information of the reconstructed events were written to Root
files.

At ILC spurious events are expected to be present in each beam-crossing. At
the ILC with $\sqrt{s}=500$ GeV an average of 1.5 low $\it{P_{T}}$ hadron
events from $\gamma\gamma$ interactions is expected per bunch crossing. A
number ($\mathcal{O}(10)$ ) of electron-positron pairs from beam-beam
interactions is also expected to reach the tracking system of the detector in
each bunch crossing. The $\gamma\gamma$ interactions were generated either
with Pythia 6.422 [25] (if $M_{\gamma\gamma}$ > 2 GeV) or a dedicated
generator [26] for $\gamma\gamma$ interactions (otherwise). The electron-
positron pairs from beam-beam interactions were generated with GuineaPig. A
pool of such events were created, and random events picked from the pool were
added to each physics event during full simulation. The SGV simulation did not
implement this mechanism, instead reconstructed objects coming from overlay
were extracted from random background events and overlaid on the signal events
at ntuple level.

Final significance values were computed adding the contribution of both
polarisations. The signal production cross-section depends on the beam
polarisation. Even more so, the background levels are expected to be quite
different for the two cases, mainly because the strongly polarisation
dependent $e^{+}e^{-}\rightarrow W^{+}W^{-}$ process is one of the most
important sources of background. Because of this the samples are weighted by
the likelihood ratio (LR) test statistic [27, 28]. As the experiment is a pure
event-counting one, the LR statistic has the simple form

$N_{\sigma}=\frac{\sum^{n_{samp}}_{i=1}s_{i}\ln{(1+s_{i}/b_{i})}}{\sqrt{\sum^{n_{samp}}_{i=1}n_{i}\left[\ln{(1+s_{i}/b_{i})}\right]^{2}}}$
(1)

where $s_{i}$ and $b_{i}$ are expected number of signal and background events
in each sample, and $N_{\sigma}$ is the significance of the test, expressed as
the equivalent number of standard deviations in the Gaussian limit. Depending
on which hypothesis is tested (exclusion or discovery), $n_{i}$ is either
$s_{i}+b_{i}$ (exclusion), or $b_{i}$ (discovery).

(a) (b) (c)

Figure 2: (a) $\widetilde{\tau}$ pair production and decay to $\tau$ and
$\widetilde{\chi}^{0}_{1}$; (b) hadronic $\tau$ decay; (c) leptonic $\tau$
decay. Blue lines indicates the $\tau$, red lines indicate undetectable final
state particles, while thick black lines indicate the detectable final state
particles.

Well identifying $\tau$:s is obviously a key requirement for this analysis. To
do this, we use the DELPHI $\tau$ finder [10]. This algorithm was particularly
developed to identify $\tau$:s in low multiplicity events. It iterative builds
$\tau$ candidates from all possible combinations of charged tracks, starting
with single tracks, then in each iteration adding tracks to existing
combinations, always retaining the combination yielding the lowest mass, and
requiring that the mass is below 2 GeV. Characteristics of $\tau$ decays are
also taken into account, so that e.g. an identified muon is not allowed to be
grouped with another track. The charged track grouping is terminated when no
more groups with mass below 2 GeV can be made. Neutrals are then added to the
found candidates, as long as the mass remains below 2 GeV - any neutrals that
can not be grouped are left as belonging to the “rest-of-event” class. To
allow for tracks not being found, it is not, at this stage, required that the
groups have charge $\pm$ 1.

##  6 Signal and background

###  6.1 Signal characterisation

Assuming R-parity conservation and assuming that the $\widetilde{\tau}$ is the
NLSP, $\widetilde{\tau}$’s will be produced in pairs via $Z^{0}$/$\gamma$
exchange in the s-channel and they will decay to a ${\tau}$ and an LSP
(assuming mass differences above the mass of the ${\tau}$, as is done in this
study). The LSP, as already mentioned, is stable and weakly interacting, hence
it will leave the detector without being detected. The ${\tau}$, with a proper
lifetime of 2.9 x 10-13 s, will decay before leaving any signal in the
detectors. The only detectable activity in the signal events is therefore the
decay products of the two ${\tau}$’s. Figure 2(a) shows the diagram of the
$\widetilde{\tau}$ production and decay, Figures 2(b) and (c) show the
diagrams of the subsequent $\tau$ decays, in the hadronic and leptonic modes,
respectively. Signal events are therefore characterised by a large missing
energy and momentum, not only due to the invisible LSPs but also to the
neutrinos from both ${\tau}$-decays. Since the $\widetilde{\tau}$’s are
scalars and hence isotropically produced, these events have a large fraction
of the detected activity in the central region of the detector.

(a) (b) (c)

Figure 3: Main irreducible background sources: (a) s-channel $WW$ production
with both $W$:s decaying to $\tau$ and $\nu$; (b) t-channel $WW$ production
with both $W$:s decaying to $\tau$ and $\nu$; (c) t-channel $ZZ$ production
with one $Z$ decaying to $\tau$:s, the other to $\nu$:s. The same colour-
coding is as in Figure 2 is used, and the various $\tau$ decays modes (cf.
Figures 2(b) and (c)) are collectively represented by the hatched circles.

The $\widetilde{\tau}$’s must also be rather heavy, so they will not have a
large boost in the lab-frame, and since the LSP is also quite heavy, the
direction of the $\widetilde{\tau}$ does not strongly correlate to that of the
visible $\tau$ after the decay. As a consequence the two $\tau$-leptons are
expected to go in directions quite independent of each other resulting in
events with un-balanced transverse momentum, large angles between the two
$\tau$-lepton directions and zero forward-backward asymmetry. These properties
are however not necessarily present in any event - the two $\tau$’s could
accidentally happen to be back to back, for example.

###  6.2 Main background sources

The main sources of background, given the generic signal topology, i.e. two
$\tau$’s and an unseen recoil system, are SM processes with real or fake
missing energy. They can be classified into “irreducible” and “almost
irreducible” sources. The first are events with two $\tau$’s and real missing
energy, i.e. neutrinos. The main contribution to this group are $ZZ$ events
with one $Z$ decaying to two neutrinos and the other to two $\tau$’s, and
fully leptonic $WW$ events, where both the $W$’s decays to $\tau$ and
neutrino. The dominating diagrams are shown in Fig. 3. In principle, $ZWW$ and
$ZZZ$ events decaying to two $\tau$’s and four neutrinos also contribute, but
are of less concern due to their low cross sections. Nevertheless, these
process are included in the simulated background samples.

(a) (b)

Figure 4: Main background sources: (a) $\tau$ pair production and decay; (b)
$\tau$ pair production in $\gamma\gamma$ events. In (b), the thick dashed red
lines indicate that while the outgoing electrons or positrons are in principle
detectable, this diagram contributes to the background only in the case that
they are deflected so little that they escape detection by remaining inside
the beam-pipe, and hence contribute to (fake) missing energy.

The second group of events are those which are not really two $\tau$’s and
neutrinos, but after reconstruction look very similar. They are events with
two soft $\tau$-jets, with two other leptons plus true missing energy or two
$\tau$’s plus fake missing energy. The main sources for events with true
missing energy in this group are on one hand $\tau$ pair production, with the
$\tau$’s decaying such that most energy goes to the neutrinos, (Figure 4(a)),
and on the other hand, $ZZ$ or $WW$ decaying to two neutrinos and at least one
lepton other than $\tau$’s, i.e. the type of processes shown in Figure 3, but
with one, or both, $\tau$’s replaced by muons or electrons. The background
with fake missing energy comes mainly from $\tau$ pair production with Initial
State Radiation (ISR) at very low angles, events with two $\tau$’s and two
very low angle electrons (below the acceptance of the detector) in the final
state (Figure 4(b)) and events where two $\tau$’s are produced by a
$\gamma\gamma$ interaction and not from an $e^{+}e^{-}$ one ; in the latter
case there is not really missing energy but the assumption that the initial
energy is the energy of the incoming electron and positron is wrong.

##  7 Analysis

###  7.1 General cuts

Taking into account the signal signature and the main background sources,
different cuts have been designed in order to separate the signal from the
background.

(a) (b)

Figure 5: Distributions for background and signal (the latter on an arbitrary
scale). (a): the number of charged particles. For clarity, the $\gamma\gamma$
background has been suppressed in this figure. (b) :
$q_{jet}\cos{\theta_{jet}}$ after the cut on the number of charged particles.
Only the irreducible $WW\rightarrow l\nu l\nu$ is shown, since this background
is what this cut is designed to reduce. Events above the vertical lines are
accepted.

The study was focused on small differences between the $\widetilde{\tau}$ and
LSP masses666Larger mass differences were also analysed in order to cross-
check and try to improve the limits from the previous studies., $M_{\tau}$ $<$
$\Delta M$ $<$ 11 GeV. This class of signal events are quite similar to high
cross-section multi-peripheral $\gamma\gamma$ events, with one important
difference, namely that in the $\gamma\gamma$ events the final state also
includes the beam-remnant electrons and positrons. These are scattered by some
low angle, so that demanding absence of signal in the calorimeter closest to
the beam pipe (the BeamCal) strongly reduces this source of background, and
was required as a pre-selection step before applying the following cuts. To
reject tracks from low invariant mass $\gamma\gamma$ events and from beam-beam
interactions (“overlay tracks”) occurring in the same bunch-crossing as the
physics event, a selection of tracks is performed, details of which are given
in Sect. 7.4.

All quantities mentioned in the following are to be understood as being
calculated using reconstructed objects passing these cuts. Tables 2 and 3 show
the cut-flow for the cuts applied at two model-points in the order the cuts
are applied in the discussion that follows. All figures in this section will
show background and signal distributions for various quantities that cuts are
applied to at the point where the cut on the shown quantity will be the next
one to be done, i.e. for events passing all previous cuts.

Initially, cuts are applied on properties that the signal often does not have,
while some backgrounds do have: The multiplicity of the event can be
constrained taking into account that the visible part of the
$\widetilde{\tau}$ decays comes only from the decays of the two $\tau$’s and
maybe an ISR photon. For that reason the number of charged particles in the
event is required to be between 2 and 10. This cut removes practically all
hadronic background, see Fig. 5(a). $WW$ events with each of the $W$’s
decaying to a lepton and a neutrino are highly forward-backward asymmetric;
they can be effectively removed by requiring the sum of the product of the
charge and the cosine of the polar angle, $q\cos{\theta_{jet}}$, of the two
most energetic jets to be above -1.0, see Fig. 5(b). $ZZ$ events with one $Z$
decaying to two neutrinos and the second one to a electron or muon pair are
highly suppressed demanding $|M_{vis}-M_{Z}|>4$ GeV, see Fig. 6 (a).

The cuts mentioned in the previous paragraph does not depend on the properties
of the signal model point. However, all following cuts applied will depend on
the model-point considered, and so will depend on either, or both, of
$M_{\widetilde{\tau}}$ and $M_{LSP}$.

(a) (b)

Figure 6: Distributions of the visible mass (a) and the missing energy (b),
for two signals (arbitrary scale), and background, after previous cuts have
been applied (cf. Tabs. 2 and 3). Events below (above) the vertical lines are
accepted in a(b) at the model points indicated by the colour of the lines. In
the $M_{vis}$ distribution, also events between the dashed lines around
$M_{Z}$ were rejected.

The first group of cuts contain those in properties that the
$\widetilde{\tau}$-events must have. Since the two LSPs from the
$\widetilde{\tau}$-decays are invisible to the detector, signal events have to
have a missing energy, $E_{miss}$, greater than $2\times M_{LSP}$.

(a) (b)

Figure 7:  Distributions of the missing transverse momentum for (a)
M${}_{\widetilde{\tau}}$ = 230 GeV $\Delta(M)$ = 34 GeV, and (b)
M${}_{\widetilde{\tau}}$ = 245 GeV $\Delta(M)$ = 10 GeV. The two signals are
on arbitrary scale, and all previous cuts have been applied (cf. Tabs. 2 and
3). Events above the vertical lines are accepted.

Likewise, the visible mass, $M_{vis}$ can not be greater than $E_{CMS}-2\times
M_{LSP}$. Therefore, events should fulfil $M_{vis}<E_{CMS}-2\times M_{LSP}$
(Fig. 6 (a)) and $E_{miss}>2\times M_{LSP}$ (Fig. 6 (b)) to be considered
further. Also a cut in the maximum total momentum, smaller than 70$\%$ the
beam momentum, is applied for the same reason: $P_{tot}<0.70E_{beam}$. There
should only be 2 or 3 clusters found by the Jade algorithm (with
$y_{cut}$=0.02) and a total seen event charge should be between -1 and 1. A
specific algorithm for $\tau$-identification was also applied. This algorithm
consists in a first set of conditions requiring to have a pattern of charged
tracks typical for $\tau$-decay, viz. exactly two $\tau$-jets (obtained with
the DELPHI tau-finder) with 1 or 3 charged particles in each charged jet, jet-
charge $\pm$1, and opposite charge between both jets. These two jets could be
accompanied by further particles, that are not compatible with the
requirements imposed by the DELPHI tau-finder to be considered as $\tau$-jets.
A further set of conditions on the jets is devoted to the reduction of
background from sources with leptons not from $\tau$-decays. The background of
single $W$, with the $W$ decaying to $\tau$ and neutrino, see Figure 8, can be
reduced by a cut depending on beam-polarisation: this process can only occur
if at least one of the beams has the correct polarisation. Since the degree of
polarisation of the positron beam is lower than that of the electron beam,
this background is more likely to yield an electron as the beam-remnant for
$\mathcal{P}_{+-}$, a positron for $\mathcal{P}_{-+}$. We therefore remove
events with an electron(positron) $\tau$ candidate in the $\mathcal{P}_{+-}$ (
$\mathcal{P}_{-+}$) samples.

Figure 8:  Single $W$ production with the $W$ decaying to $\tau$ and $\nu$.
The labels indicate the only beam-polarisations for which this process is
possible: either left-handed electrons or right-handed positrons.

This background together with the background from $WW\rightarrow
e\nu_{e}\mu\nu_{\mu}$ and from $\gamma\gamma$ events with a beam-remnant
deflected to larger angles is further reduced by rejecting those events in
which the most energetic jet consists of a single electron. The two charged
jets were also required to not be made by single leptons with the same
flavour. These selections reduce the signal efficiency to 38$\%$ but with a
reduction of the background of the order of 94$\%$, depending on the region of
the SUSY parameter space.

Since the $\widetilde{\tau}$-decay is a two body decay, it is possible to
determine the maximum and minimum momentum of each of the decay products as a
function of the $\widetilde{\tau}$ mass, the mass of LSP and the centre-of-
mass energy of the collider. A cut in the minimum momentum can not be applied
due to the presence of neutrinos in the $\tau$ decay, with the corresponding
decrease of observable momentum. The maximum value can be used even if it is
smeared by the missed neutrinos. The expression for the maximum jet momentum
is given by:

$P_{max}=\frac{\sqrt{s}}{4}\left(1-\left(\frac{M_{LSP}}{M_{\widetilde{\tau}}}\right)^{2}\right)\left(1+\sqrt{1-\frac{4M_{\widetilde{\tau}}\leavevmode\nobreak\
^{2}}{s}}\right)$ (2)

and events with the higher jet momentum above this limit are excluded from the
further analysis at the given model point. Excluding the cut applied by the
$\tau$-identification algorithm, the signal efficiency for each of the cuts is
at least 95$\%$ at all model points.

(a) (b)

Figure 9:  Distributions of the acoplanarity angle between the two jets for
(a) M${}_{\widetilde{\tau}}$ = 230 GeV $\Delta(M)$ = 34 GeV, and (b)
M${}_{\widetilde{\tau}}$ = 245 GeV $\Delta(M)$ = 10 GeV. The two signals are
on arbitrary scale, and all previous cuts have been applied (cf. Tabs. 2 and
3). Events between the vertical lines are accepted.

A second group of cuts is based on those properties that the
$\widetilde{\tau}$-events might have, but will rarely be

(a) (b)

Figure 10:  Distributions of the polar angle of the vectorial sum of the two
jet momenta for (a) M${}_{\widetilde{\tau}}$ = 230 GeV $\Delta(M)$ = 34 GeV,
and (b) M${}_{\widetilde{\tau}}$ = 245 GeV $\Delta(M)$ = 10 GeV. The two
signals are on arbitrary scale, and all previous cuts have been applied (cf.
Tabs. 2 and 3). Events below the vertical lines are accepted.

present in background events. As already pointed out, the $\widetilde{\tau}$’s
are scalars, and therefore isotropically produced, while the backgrounds are
either fermions or vector bosons, and tend to be produced at small angles to
the beam axis. This allows to impose cuts requiring events to have high
missing transverse momentum ($P_{Tmiss}$), and large acoplanarity
$\theta_{acop}$. The distributions for these two quantities are shown in Figs.
7 and 9, for two model-points as they are after applying all previous cuts,
applicable to the model-points in question. The total seen momentum,
$\sum\bar{p}$, tends to be in an almost random direction in signal events,
while it tends to point close to the beam-axis for most backgrounds. Therefore
a cut on the direction of the total momentum,
$|\cos(\theta_{\sum\bar{p}})|<0.85$, was imposed, see Fig.10. In addition, a
cut on the variable $\rho$ is imposed. This variable is calculated by first
projecting the event on the x-y (transverse) plane, and calculating the thrust
axis in that plane. $\rho$ is then the transverse momentum (in the plane) with
respect to the thrust axis. The cut in $\rho$ helps to reject events with two
$\tau$’s back-to-back in the transverse projection with the visible part of
the decay of one of the $\tau$’s in the direction of its parent, while the
other $\tau$ decays with the invisible $\nu$ closely aligned with the
direction of its parent. These events fake the signal topology, having both a
large missing transverse momentum and high acoplanarity. However, they would
have a small value of $\rho$. The distributions of $\rho$ for signal and
background, at this point in the cut chain, are shown in Fig. 11, for two
model points. The values at which the cuts described in this paragraph are set
depends on the mass difference between the $\widetilde{\tau}$ and the LSP
(except the cut on $\theta_{\sum\bar{p}}$), and are detailed in Tab. 1.

Table 1: Cuts depending on the mass difference between the $\widetilde{\tau}$ and the LSP. The specific cuts applied for $\Delta(M)$ = 2 GeV are explained in Sect. 7.2. $\Delta(M)$ | $P_{Tmiss}$ | $\rho$ | $\theta_{acop}$
---|---|---|---
[GeV] | [GeV] | [GeV] | [Rad]
34 |  $>$ 11 |  $>$ 11 | $\in[0.2,3.1]$
10 |  $>$ 4 |  $>$ 8 | $\in[0.2,2.0]$
3 |  $>$ 2 |  $>$ 3 | $\in[0.2,2.3]$
2 | - | no cut | $\in[0.2,2.3]$

(a) (b)

Figure 11:  Distributions of the variable $\rho$, described in the text, for
(a) M${}_{\widetilde{\tau}}$ = 230 GeV $\Delta(M)$ = 34 GeV, and (b)
M${}_{\widetilde{\tau}}$ = 245 GeV $\Delta(M)$ = 10 GeV. The two signals are
on arbitrary scale, and all previous cuts have been applied (cf. Tabs. 2 and
3). Events above the vertical lines are accepted.

Cutting in these properties has a certain cost in efficiency but improves the
signal-to-background ratio.

As a final cut it is noted that sizeable energy detected at low angles to the
beam is rare in signal events, but common in many backgrounds. Events with
more than 2 GeV detected at angles lower than 20 degrees to the beam axis are
therefore rejected. This cut is however not useful for small mass differences.
After applying these cuts the main sources of remaining background are $WW$
events with each $W$ decaying to $\tau\nu$ and events with four fermions in
the final state coming from $\gamma\gamma$ interactions, mostly $\tau\tau$
events. As an illustration of the expected signal after all cuts, Fig 12 shows
the distribution of the higher of the two jet-momenta
($p_{jet\leavevmode\nobreak\ max}$) for the model-point with
M${}_{\widetilde{\tau}}$ = 230 GeV and $\Delta(M)$ = 34 GeV, when all
described cuts except the one on $p_{jet\leavevmode\nobreak\ max}$ has been
applied. In this figure, the signal and background distributions correspond to
the same 1600 fb-1 of integrated luminosity.

Table 2: Flow table M${}_{\widetilde{\tau}}$ = 230 GeV $\Delta(M)$ = 34 GeV. The numbers correspond to an integrated luminosity of 1.6 ab-1 and polarisation $\mathcal{P}_{+-}$. Cut | Signal | $e^{+}e^{-}\rightarrow$ | $e^{+}e^{-}\rightarrow$ | $e^{+}e^{-}\rightarrow$ | $\gamma\gamma(\gamma e)\rightarrow$ | All other
---|---|---|---|---|---|---
|  | $\tau\tau/\mu\mu\leavevmode\nobreak\ \nu\nu$ | $ll\leavevmode\nobreak\ \nu\nu$ | $\tau\tau/\mu\mu$ | 2(3) fermions |
|  |  | ($\geq 1\leavevmode\nobreak\ l$ not $\tau$) | (no $\nu$’s) |  |
no cut | $8292.8$ | $134\,100$ | $533\,900$ | 4.377 $\cdot 10^{6}$ | 1234 $\cdot 10^{6}$ | 2616 $\cdot 10^{6}$
BeamCal veto | $8284.5$ | $133\,966$ | $532\,800$ | 4.368 $\cdot 10^{6}$ | 1086 $\cdot 10^{6}$ | 2603 $\cdot 10^{6}$
2 $\leq N_{charged}\leq$ 10 | $8219.0$ | $108\,900$ | $403\,100$ | 3.542 $\cdot 10^{6}$ | 1024 $\cdot 10^{6}$ | 1720 $\cdot 10^{6}$
$q\cos{\theta_{jet}}>-1$ | $7179.1$ | $70\,720$ | $216\,800$ | 2.671 $\cdot 10^{6}$ | 841.6 $\cdot 10^{6}$ | 1504 $\cdot 10^{6}$
$|M_{vis}-M_{Z}|>4$ | $7178.2$ | $66\,820$ | $207\,200$ | 2.552 $\cdot 10^{6}$ | 826.8 $\cdot 10^{6}$ | 1500 $\cdot 10^{6}$
$E_{miss}>2M_{LSP}$ | $7177.4$ | $19\,620$ | $14\,570$ | $365\,493.0$ | 238.0 $\cdot 10^{6}$ | 1399 $\cdot 10^{6}$
$M_{vis}<\sqrt{s}-2M_{LSP}$ | $7103.6$ | $17\,990$ | $12\,490$ | $334\,886.7$ | 235.5 $\cdot 10^{6}$ | 1398 $\cdot 10^{6}$
$\tau$ Id | $2953.9$ | $1250$ | $451.2$ | $29\,088.3$ | 23.96 $\cdot 10^{6}$ | 166.9$\cdot 10^{6}$
$P_{jet}<P_{max}$ | $2953.9$ | $1003$ | $268.3$ | $23\,963.2$ | 23.55 $\cdot 10^{6}$ | 165.6$\cdot 10^{6}$
$P_{Tmiss}$ > 11 GeV | $2478.7$ | $242.5$ | $175.6$ | $3481$ | $83\,970$ | $54.7$
$\theta_{acop}\in$ [0.2, 3.1] | $2293.8$ | $208.6$ | $144.9$ | $3075$ | $69\,190$ | $52.0$
$|\cos(\theta_{\sum\bar{p}})|<0.85$ | $2125.4$ | $147.1$ | $113.6$ | $1876$ | $32\,520$ | $39.8$
$\rho$ > 11 | $1317.7$ | $60.3$ | $39.3$ | $0.0$ | $1923$ | $19.3$
$E_{<20^{\circ}}<2$ GeV | $1015.9$ | $39.6$ | $27.2$ | $0.0$ | $81.0$ | $12.1$
Table 3: Flow table M${}_{\widetilde{\tau}}$ = 245 GeV $\Delta(M)$ = 10 GeV. The numbers correspond to an integrated luminosity of 1.6 ab-1 and polarisation $\mathcal{P}_{+-}$. Cut | Signal | $e^{+}e^{-}\rightarrow$ | $e^{+}e^{-}\rightarrow$ | $e^{+}e^{-}\rightarrow$ | $\gamma\gamma(\gamma e)\rightarrow$ | All other
---|---|---|---|---|---|---
|  | $\tau\tau/\mu\mu\leavevmode\nobreak\ \nu\nu$ | $ll\leavevmode\nobreak\ \nu\nu$ | $\tau\tau/\mu\mu$ | 2(3) fermions |
|  |  | ($\geq 1\leavevmode\nobreak\ l$ not $\tau$) | (no $\nu$’s) |  |
no cut | $672.8$ | $134\,100$ | $533\,900$ | 4.377 $\cdot 10^{6}$ | 1234 $\cdot 10^{6}$ | 2616 $\cdot 10^{6}$
BeamCal veto | $672.1$ | $133\,966$ | $532\,800$ | 4.368 $\cdot 10^{6}$ | 1086 $\cdot 10^{6}$ | 2603 $\cdot 10^{6}$
2 $\leq N_{charged}\leq$ 10 | $664.3$ | $106\,400$ | $403\,100$ | 3.542 $\cdot 10^{6}$ | 1024 $\cdot 10^{6}$ | 1720 $\cdot 10^{6}$
$q\cos{\theta_{jet}}>-1$ | $585.2$ | $70\,720$ | $216\,800$ | 2.671 $\cdot 10^{6}$ | 841.7 $\cdot 10^{6}$ | 1504 $\cdot 10^{6}$
$|M_{vis}-M_{Z}|>4$ | $585.2$ | $66\,820$ | $207\,200$ | 2.552 $\cdot 10^{6}$ | 826.8 $\cdot 10^{6}$ | 1500 $\cdot 10^{6}$
$E_{miss}>2M_{LSP}$ | $585.1$ | $6218$ | $886.6$ | $120\,800$ | 179.0 $\cdot 10^{6}$ | 1043 $\cdot 10^{6}$
$M_{vis}<\sqrt{s}-2M_{LSP}$ | $583.9$ | $5669$ | $758.0$ | $109\,100$ | 172.7 $\cdot 10^{6}$ | 969.9 $\cdot 10^{6}$
$\tau$ Id | $230.6$ | $484.0$ | $56.3$ | $13\,020$ | 1.010 $\cdot 10^{6}$ | 127.4 $\cdot 10^{6}$
$P_{jet}<P_{max}$ | $229.9$ | $444.2$ | $38.2$ | $11\,150$ | $822\,400$ | 108.0 $\cdot 10^{6}$
$P_{Tmiss}$ > 4 GeV | $177.4$ | $33.4$ | $5.1$ | $435.0$ | $36\,210$ | $84\,660$
$\theta_{acop}\in$ [0.2, 2.0] | $164.6$ | $30.0$ | $5.1$ | $299.6$ | $32\,310$ | $71\,650$
$|\cos(\theta_{\sum\bar{p}})|<0.85$ | $156.4$ | $22.7$ | $0.0$ | $290.0$ | $24\,220$ | $20\,290$
$\rho$ > 8 | $55.6$ | $2.7$ | $0.0$ | $0.0$ | $798.6$ | $0.5$
$E_{<20^{\circ}}<2$ GeV | $55.6$ | $2.7$ | $0.0$ | $0.0$ | $610.2$ | $0.5$

The polarisation plays an important role in the capability of
excluding/discovering the different regions of the SUSY space. Table 4 shows
the number of signal and background events for a specific model point for the
two main ILC running polarisations and for unpolarised beams. The difference
in the number of signal events arises from the dependence of the cross section
on the polarisation, as well as from the effect on selection efficiency due to
the polarisation of the $\tau$ coming from the $\widetilde{\tau}$. The
dependence of the cross section on the polarisation is the main factor for the
difference in $WW$ events, $e^{+}e^{-}\rightarrow\tau\nu\tau\nu$. One can see
that the signal-to-background ratio is clearly enhanced for the case of mainly
right-handed electrons, left-handed positrons. Taking the definition of
exclusion at 95$\%$ CL as $N_{\sigma}>2$ (cf. Eq. 1)777Note that Eq. 1
trivially states that $N_{\sigma}=S/\sqrt{S+B}$ in the case that there is only
one sample., it can be seen that combining the two polarisation samples with
the LR statistic does strengthen the limit, even though the $\mathcal{P}_{-+}$
sample gives a much weaker limit than the $\mathcal{P}_{+-}$ one. It is also
shown that unpolarised beams would allow neither exclusion nor discovery, even
though the sample corresponds to the same integrated luminosity as the sum of
the two polarised samples. This is not only because of the possibility to
combine samples with different beam polarisations in an optimal way is not
available for unpolarised beams; it is also because the ILC has both beams
polarised, meaning that more than half of the collisions are between opposite
polarised particles, which is necessary to allow for s-channel processes, such
as $\widetilde{\tau}$ pair production. Polarisation is not only important in
the enhancement of the signal over background but also plays an important role
in the parameter determination.

Table 4: Remaining signal and background events after the application of the selection cuts for $M_{\widetilde{\tau}}$=245 GeV and mass difference with the LSP of 10 GeV.  Polarisation | Signal | $e^{+}e^{-}\rightarrow$ | $\gamma\gamma(\gamma e)\rightarrow$ | All other | Nσ
---|---|---|---|---|---
|  | $\tau\tau/\mu\mu\leavevmode\nobreak\ \nu\nu$ | 2(3) fermions |  |
$\mathcal{P}_{+-}$ | $55.6$ | $2.7$ | $610.2$ | $0.5$ | $2.1$
$\mathcal{P}_{-+}$ | $36.0$ | $39.6$ | $1208$ | $0.5$ | $1.0$
Combined | - | - | - | - | $2.4$
Unpolarised | $73.9$ | $34.1$ | $1818$ | $1.0$ | $1.7$
(a) Figure 12:  Distribution of the higher of the two jet momenta after all
cuts (except the one on $p_{jet\leavevmode\nobreak\ max}$ itself)
corresponding to the model point M${}_{\widetilde{\tau}}$ = 230 GeV,
$\Delta(M)$ = 34 GeV. The maximum $P_{jet}$ possible at this model point is
47.63 GeV (Eq. 2), and is indeed the endpoint of the signal distribution.
(Contrary to other figures in this section, the signal is stacked on top of
the background, and corresponds to the same integrated luminosity as the
background.)

###  7.2  Low $\Delta(M)$ cuts

The cuts described above are suited for mass differences down to 3 GeV. When
the mass difference is between 3 GeV and the mass of the $\tau$ 888For mass
differences below the mass of the $\tau$ the lifetime of the
$\widetilde{\tau}$ increases exponentially and the study has to be done based
on a signature of long-lived particles travelling through the detector. the
kinematics of the signal events is very close to that of the $\gamma\gamma$
background events and the described cuts are not enough for
discovering/excluding the signal. An additional cut was done based on the
Initial State Radiation photons (ISR). Events with isolated photons with
sizeable energy were selected, allowing to extension of the limits into the
region under study. To select these, the photon should be at above 7∘ to the
beam - the lower edge of the acceptance of the tracking system of ILD - to
avoid that electrons would be mistaken for photons, and it should have an
energy above 1.1 GeV.

This cut is effective against the remaining $\gamma\gamma$ background because
these events become candidates due to fake missing transverse momentum. If the
presence of an ISR is requested, the incoming electron or positron that
emitted the ISR must have recoiled against the ISR. Since this is a scattering
process, not an annihilation one, the electron (positron) is still present in
the final state. Therefore, if it is required to see a high transverse
momentum ISR, the final state electron (positron) will have acquired a recoil
transverse momentum big enough to be deflected into the BeamCal, and thus to
have been rejected already at the pre-selection stage. On the other hand, if
the ISR was emitted from an electron or positron that was subsequently
annihilated into a $Z$, as is the case for the signal process, the transverse
momentum of the ISR is included in the decay products of the $Z$, and no
signal is expected in the BeamCal.

Events containing such ISR photons, and passing all cuts except those detailed
in Tab. 1, were retained if $P_{Tmiss}>4$ GeV and the sum of the absolute
value of the jet momenta was between 1 and 4 GeV.

###  7.3 Determining the worst scenario

The aim of this study is to evaluate the capabilities of the
$\widetilde{\tau}$ search at ILC with no loopholes. One of the properties of
the $\widetilde{\tau}$ \- apart from its and the LSP masses - is the
$\widetilde{\tau}$ mixing. This enters in two ways: The production cross-
section depends on it, as does the polarisation of the $\tau$. The dependence
of the cross-section on the $\widetilde{\tau}$ mixing angle is shown in Fig.
13, for the two beam-polarisations considered at ILC. Also shown is dependence
for unpolarised beams, which shows a minimum at 53∘. The polarisation of the
$\tau$ is also important, since it influences the momentum distribution of the
${\tau}$-decay products, and hence the signal efficiency: if the $\tau$ is
pre-dominantly left-handed, and since there are no right-handed neutrinos in
the standard model, the invisible neutrino will tend to align with the
direction of the $\tau$ and take a larger fraction of the momentum, and
consequently the visible system will take less. This is illustrated in Fig. 14
which shows the momentum distribution of the pions coming from ${\tau}$-decays
(in the single pion decay mode) for different $\widetilde{\tau}$ mixing angles
and a bino-LSP. The ${\tau}$ polarisation from the $\widetilde{\tau}$ decay
depends not only on the $\widetilde{\tau}$ but also on the neutralino mixing
angle. This is because of the different behaviour of the interaction between
higgsinos and binos with the $\widetilde{\tau}$. The higgsino (like the Higgs)
carries weak hyper-charge, while the bino (like the $Z$ or the photon) does
not. Therefore, with a higgsino LSP, the $\widetilde{\tau}$ produces a
${\tau}$ with the opposite chirality with respect to that of the
$\widetilde{\tau}$, while with a bino LSP, the produced ${\tau}$ has the same
chirality. Since the lowest efficiency is expected for left-handed $\tau$’s,
we evaluate the efficiency for the assumption that the neutralino is a pure
Bino for $\widetilde{\tau}$ mixings below 45∘ (i.e. for a $\widetilde{\tau}$
more left than right), and for a pure higgsino LSP for mixings above 45∘. This
assures that the most difficult case is studied at any mixing. In Fig. 15(a)
the signal efficiency for an example point is shown, with these assumptions.

Figure 13:  Cross section for $\widetilde{\tau}$ pair production as a function
of the $\widetilde{\tau}$ mixing angle for the two main ILC polarisations. The
green dotted line shows the case corresponding to the simply combination of
both samples, corresponding to unpolarised beams for equal integrated
luminosity.  Figure 14: Momentum distribution of the pions coming from
${\tau}$-decays in the mode $\tau\rightarrow\pi^{\pm}\nu$ for different
$\widetilde{\tau}$ mixing angles. The neutralino was assumed to be a pure
bino.

(a) (b)

Figure 15: (a): Signal efficiency as a function of the $\widetilde{\tau}$
mixing angle for both main ILC polarisations. The LSP composition was chosen
giving the worst case for each mixing. (b):Number of signal and background
selected events as a function of the $\widetilde{\tau}$ mixing angle for both
main ILC polarisations.

Combining cross sections and signal efficiencies, the number of selected
signal events for each polarisation as a function of the $\widetilde{\tau}$
mixing angle was evaluated. Figure 15(b) shows the obtained values together
with the number of background events passing the selection cuts, which is
different for each beam polarisation but obviously do not have any dependence
on the $\widetilde{\tau}$ mixing angle. As can be seen, the background is much
larger for the $\mathcal{P}_{-+}$ sample. This is because the irreducible
background from $e^{+}e^{-}\rightarrow WW\rightarrow\tau\nu\tau\nu$ is
strongly polarisation dependent. This constitutes the final ingredient to
estimate the worst-case scenario: from these values, the signal over
background significance as a function of the mixing angle was computed for
each polarisation, as shown in figure 16(a). Final significance values were
computed adding the contribution of both polarisations weighted by the
likelihood ratio statistic. The results are plotted in figure 16(b). One can
see that rather uniform sensitivity to all mixing angles is obtained, and that
for the smallest mass differences - the ones closest to the critical region -
a mixing angle around 53∘ can indeed be considered as the worst one, and
validates our choice of this mixing angle being the worst case999 Note that
the estimation of the worst scenario was not done with the final cuts used in
the analysis. The final cuts were further optimised taking into account the
whole parameter space. This does not affect the conclusion of this section..

(a) (b)

Figure 16: (a): Signal significance as a function of the $\widetilde{\tau}$
mixing angle for both main ILC polarisations. (b)Signal significance weighting
both polarisations using the likelihood-ratio static in the H20 ILC
conditions.

###  7.4 Effect of overlay tracks

The overlay tracks can not be neglected in the $\widetilde{\tau}$ studies,
since they have similar properties to the visible tracks from the
$\widetilde{\tau}$-decays in the region with small mass differences between
$\widetilde{\tau}$ and LSP. The main characteristics of the overlay tracks are
the low transverse momentum and low angles to the beam axis. The tracks in low
$\it{P_{T}}$ hadron events usually originates in the beam-spot, contrary to
the tracks from the $\tau$ decays. Figure 17(a) compares transverse momentum
distributions for signal and overlay tracks, and it can be seen that while the
overlay tracks can easily be removed by a cut in this quantity for large to
moderate mass-differences, this is not possibly for the smallest ones. Figure
17(b) makes the same comparison for the cosine of the polar angle of the
tracks, showing that in this variable, the difference between overlay and
signal tracks does not depend on the mass-difference.

(a) (b)

Figure 17: (a): Transverse momentum distribution for signal and overlay
tracks. Signal tracks correspond to different mass differences between
$\widetilde{\tau}$ and LSP. (b): Distribution of cosine of polar angle for
signal and overlay tracks. Signal tracks correspond to different mass
differences between $\widetilde{\tau}$ and LSP.

In the type of event-topologies the signal gives rise to, viz. two displaced
$\tau$ decays, no primary vertex can be fitted. However, at ILC the transverse
beam-spot size is significantly smaller than the measurement precision, so
that in that plane, the primary vertex position is known, even without
fitting. The longitudinal size of the beam-spot, on the other hand, is much
larger than the measurement precision. Figure 18 shows the significance of the
impact parameter with respect to the origin in the transverse direction for
overlay and signal tracks for different mass differences. Here, the impact
parameter has been life-time signed, i.e. the sign is determined by the
position w.r.t. to the beam-spot where the track intersects the axis of the
jet it belongs to, taking the direction of the jet as the positive direction.
For the overlay tracks, the expected insignificant transverse impact
parameters for tracks originating inside the beam-spot is observed. For the
signal case, one observes more and more significant deviation in the impact
parameters, reflecting that the tracks do not originate from inside the beam-
spot.

(a) overlay (b) $\Delta(M)$ = 2 GeV (c) $\Delta(M)$ = 10 GeV (d) $\Delta(M)$ =
34 GeV

Figure 18: Transverse impact parameter significance of tracks from overlay and
from signals with $M_{\widetilde{\tau}}$ = 240 GeV, and different mass-
differences.

Hence, for model points with mass differences $\leq$ 10 GeV, we select tracks
with $|cos{\theta}|<0.7$, and the transverse impact parameter significance
larger than 2. If the mass difference is $>$ 10 GeV, we instead require that.
$P_{T}>2$ GeV.

While these cuts clearly occasionally will remove tracks from the $\tau$ decay
resulting in the event being lost as a candidate (since the strict topology
requirements of exactly two $\tau$-jets with the correct charges will not be
met), this is compensated by previous lost candidates being salvaged by
removing wrongly included tracks from overlay in the topology evaluation.
This, together with the fact that many background events only became
candidates because of the presence of overlay tracks, and so will be discarded
by the cuts, results in a net amelioration of the sensitivity. Figure 19 shows
the results obtained with and without cuts together with the results obtained
if the overlay tracks were added neither to signal nor background. For the
case with the smallest mass difference, shown in Figure 19(a), there is a
strong reduction of the significance when adding overlay tracks. For the
larger mass-difference (Figure 19(b)), there is only a slight degradation due
to the overlay tracks. For both mass differences, the overlay removal
procedure indeed ameliorates the sensitivity.

(a) $\Delta(M)=$ 3 GeV (b) $\Delta(M)=$ 10 GeV

Figure 19: Number of sigmas for a $\widetilde{\tau}$ with mass 240 GeV and
different mass-differences. The plot assumes H20 ILC scenario combining both
polarisations using the likelihood-ratio statistic. Blue lines correspond to
the case with all the tracks and the red ones after rejecting tracks not
satisfying the cuts described in the text. The green curves correspond to the
study without overlay tracks.

##  8 Exclusion/discovery limits

The exclusion and discovery limits extracted from this study are shown in
Figure 20. They assume the lightest $\widetilde{\tau}$ to be the NLSP and the
lightest neutralino the LSP, and are valid for any $\widetilde{\tau}$ mixing
angle. It is also relevant to compare these results with the current
$\widetilde{\tau}$ limits coming from LEP. The LEP limits are also valid for
any value of the not shown model parameters.

Figure 20: Exclusion and discovery $\widetilde{\tau}$ limits from the current
studies compared to the ones from LEP studies [8]. The ILC region corresponds
to the full H20 data-set at $E_{CMS}$ = 500 GeV, i.e. 1.6 ab-1 at each of the
beam polarisations $\mathcal{P}_{+-}$ and $\mathcal{P}_{-+}$.

The projection of the limits in the $M_{\widetilde{\tau}}$-$\Delta$M plane is
shown in figure 21. The region for mass differences below the mass of the
$\tau$, not included in the current study, is shown for completeness. In the
region with $\Delta$M larger than $M_{\tau}$ exclusion and discovery ILC
limits are compared to the ones from LEP. The projected HL-LHC limits are also
shown. Since they are highly model-dependent, the comparison in this case have
to be taken with care: here limits considering only the
$\widetilde{\tau}_{R}$-pair production are shown, since, while still being
optimistic, they are closest to the ones expected for the lightest
$\widetilde{\tau}$ at minimal cross-section. It should be noted that the HL-
LHC projection is only an exclusion limit - no discovery potential is
expected.

For the region with $\Delta$M smaller than $M_{\tau}$ results from LEP and LHC
are shown. The LEP studies cover not only the region where the
$\widetilde{\tau}_{1}$ travels through the detector without decaying but also
the region with decays at a certain distance from the production vertex. In
those regions acoplanar leptons, tracks with large impact parameters and
kinked tracks are looked for, depending on the $\widetilde{\tau}_{1}$ lifetime
[29, 30]. The figure also shows the extrapolation of the ILC limits for the
scenarios with centre-of-mass energy 250 GeV and 1 TeV.

Figure 21: $\widetilde{\tau}$ limits in the $M_{\widetilde{\tau}}$-$\Delta$M
plane. ILC results from the current studies are shown together with limits
from LEP [8] and projection for HL-LHC [9]. The ILC region corresponds to the
full H20 data-set at $E_{CMS}$ = 500 GeV, i.e. 1.6 ab-1 at each of the beam
polarisations $\mathcal{P}_{+-}$ and $\mathcal{P}_{-+}$. The region with mass
differences below the mass of the ${\tau}$ is also shown with LEP [29, 30] and
LHC results [31], even if it is not covered by this study. In addition, the
extrapolation of the ILC current results to the ILC 250 GeV and 1 TeV running
scenarios is shown.

###  8.1 Comparison with previous study

Results from previous ILC studies, which used a $\int\mathcal{L}$dt = 500 fb-1
sample of events with beam polarisation $\mathcal{P}_{+-}$ are shown in Figure
22, together with results obtained under those conditions with the present
analysis.

Figure 22:  Exclusion and discovery $\widetilde{\tau}$ limits as a function of
the $\widetilde{\tau}$ and LSP masses. Exclusion limits from previous ILC
studies [14] are also shown, as well as the current analysis applied to the
same data-sample as the previous one, i.e. only the $\mathcal{P}_{+-}$ beam
polarisation sample.

The comparison of these two curves shows that the extension of the limits is
not only due to an increase of the total integrated luminosity but also to an
improvement of the analysis. The main reason of this improvement is the
application of individual limits depending on the $\widetilde{\tau}$ mass and
the mass difference. The previous studies were only making a difference for
mass differences above or below 10 GeV and were not optimised for the low mass
difference region. Another difference in the analysis is a change in the
${\tau}$-identification algorithm, excluding events with two jets consisting
of single leptons of the same flavour.

##  9 Outlook and conclusions

The capability of the ILC for excluding/discovering $\widetilde{\tau}$-pair
production up to a few GeV below the kinematic limit, without model
dependencies and even in the worst scenario, has been shown.

The worst scenario for $\widetilde{\tau}$-pair production at the ILC was
reviewed taking into account ILC beam polarisation conditions. Equal sharing
of $P(e^{-},e^{+})=(+80\%,-30\%)$ and $P(e^{-},e^{+})=(-80\%,+30\%)$ foreseen
in H20 ensures a quite uniform sensitivity to all mixing angles.

The effect of overlay tracks on the signal/background ratio for
$\widetilde{\tau}$ searches was analysed. It was found that the effect is
considerable, both for small and moderate mass differences. Cuts to mitigate
the effect were studied and applied.

The calculation of the exclusion/discovery limits in the region with mass
differences below the ${\tau}$ mass, meaning an exponential increase of the
$\widetilde{\tau}$ lifetime and consequently a study of long-lived particles
going through or decaying in different parts of the detector, is foreseen for
a future study, as is the verification that events containing only overlay
tracks cannot fake signal.

##  10 Acknowledgements

We would like to thank the LCC generator working group for producing the Monte
Carlo samples used in this study. We also thankfully acknowledge the support
by the the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
under Germany’s Excellence Strategy EXC 2121 “Quantum Universe” 390833306.
This work has benefited from computing services provided by the German
National Analysis Facility (NAF) [32].

## References

  * [1] Stephen P. Martin  “A Supersymmetry primer”  In _Adv. Ser. Direct. High Energy Phys._ 18, 1998, pp. 1–98  DOI: 10.1142/9789812839657_0001
  * [2] J. Wess and B. Zumino  “Supergauge Transformations in Four-Dimensions” [,24(1974)]  In _Nucl. Phys._ B70, 1974, pp. 39–50  DOI: 10.1016/0550-3213(74)90355-1
  * [3] Hans Peter Nilles  “Supersymmetry, Supergravity and Particle Physics”  In _Phys. Rept._ 110, 1984, pp. 1–162  DOI: 10.1016/0370-1573(84)90008-5
  * [4] Howard E. Haber and Gordon L. Kane  “The Search for Supersymmetry: Probing Physics Beyond th e Standard Model”  In _Phys. Rept._ 117, 1985, pp. 75–263  DOI: 10.1016/0370-1573(85)90051-1
  * [5] Riccardo Barbieri, S. Ferrara and Carlos A. Savoy  “Gauge Models with Spontaneously Broken Local Supersymmetry”  In _Phys. Lett._ 119B, 1982, pp. 343  DOI: 10.1016/0370-2693(82)90685-2
  * [6] “The International Linear Collider Technical Design Report - Volume 1: Executive Summary”, 2013  arXiv:1306.6327 [physics.acc-ph]
  * [7] T. J.. and K.. Olive  “Neutralino-Stau Coannihilation and the Cosmological Upper Limit on the Mass of the Lightest Supersymmetric Particle”  In _Phys. Lett. B_ 444, 1998, pp. 367  arXiv:9810360 [hep-ph]
  * [8] LEP SUSY Working Group, ALEPH, DELPHI, L$3$ and OPAL Collaborations  “Combined LEP Selectron/Smuon/Stau Results, $183$-$208$ GeV”  URL: http://lepsusy.web.cern.ch/lepsusy/www/sleptons_summer04/slep_final.html
  * [9] ATLAS Collaboration  “Prospects for searches for staus, charginos and neutralinos at the high luminosity LHC with the ATLAS Detector”, 2018  URL: https://cds.cern.ch/record/2651927
  * [10] J. Abdallah  “Searches for supersymmetric particles in e+ e- collisions up to 208-GeV and interpretation of the resul ts within the MSSM”  In _Eur. Phys. J._ C31, 2003, pp. 421–479  DOI: 10.1140/epjc/s2003-01355-5
  * [11] Georges Aad  “Search for direct stau production in events with two hadronic $\tau$-leptons in $\sqrt{s}=13$ TeV $pp$ collisions with the ATLAS detector”  In _Phys. Rev. D_ 101.3, 2020, pp. 032009  DOI: 10.1103/PhysRevD.101.032009
  * [12] Albert M Sirunyan  “Search for direct pair production of supersymmetric partners to the $\tau$ lepton in proton-proton collisions at $\sqrt{s}=$ 13 TeV”  In _Eur. Phys. J. C_ 80.3, 2020, pp. 189  DOI: 10.1140/epjc/s10052-020-7739-7
  * [13] CMS Collaboration  “Search for supersymmetry with direct stau production at the HL-LHC with the CMS Phase-2 detector”, 2019  URL: https://cds.cern.ch/record/2647985
  * [14] Mikael Berggren  “Simplified SUSY at the ILC”  In _Community Summer Study 2013: Snowmass on the Mississippi_ , 2013  arXiv:1308.1461 [hep-ph]
  * [15] Sebastian Baum, Pearl Sandick and Patrick Stengel  “Hunting for scalar lepton partners at future electron colliders”  In _Phys. Rev. D_ 102.1, 2020, pp. 015026  DOI: 10.1103/PhysRevD.102.015026
  * [16] T. Barklow et al.  “ILC Operating Scenarios”, 2015  arXiv:1506.07830 [hep-ex]
  * [17] Halina Abramowicz  “International Large Detector: Interim Design Report”, 2020  arXiv:2003.01116 [physics.ins-det]
  * [18] Wolfgang Kilian, Thorsten Ohl and Jurgen Reuter  “WHIZARD: Simulating Multi-Particle Processes at LHC and ILC”  In _Eur. Phys. J. C_ 71, 2011, pp. 1742  DOI: 10.1140/epjc/s10052-011-1742-y
  * [19] D. Schulte  “Beam-beam simulations with GUINEA-PIG”, 1999  URL: https://cds.cern.ch/record/382453
  * [20] M. Petrič et al.  “Detector simulations with DD4hep”  In _J. Phys. Conf. Ser._ 898.4, 2017, pp. 042015  DOI: 10.1088/1742-6596/898/4/042015
  * [21] F. Gaede  “Marlin and LCCD: Software tools for the ILC”  In _Nucl. Instrum. Meth. A_ 559, 2006, pp. 177–180  DOI: 10.1016/j.nima.2005.11.138
  * [22] A. Tsaregorodtsev  “DIRAC: A community grid solution”  In _J. Phys. Conf. Ser._ 119, 2008, pp. 062048  DOI: 10.1088/1742-6596/119/6/062048
  * [23] Mikael Berggren  “SGV 3.0 - a fast detector simulation”  In _International Workshop on Future Linear Colliders (LCWS11)_ , 2012  arXiv:1203.0217 [physics.ins-det]
  * [24] Stanislaw Jadach, Johann H. Kuhn and Zbigniew Was  “TAUOLA: A Library of Monte Carlo programs to simulate decays of polarized tau leptons”  In _Comput. Phys. Commun._ 64, 1990, pp. 275–299  DOI: 10.1016/0010-4655(91)90038-M
  * [25] Torbjorn Sjostrand, Stephen Mrenna and Peter Z. Skands  “PYTHIA 6.4 Physics and Manual”  In _JHEP_ 05, 2006, pp. 026  DOI: 10.1088/1126-6708/2006/05/026
  * [26] Pisin Chen, Timothy L. Barklow and Michael E. Peskin  “Hadron production in gamma gamma collisions as a background for e+ e- linear colliders”  In _Phys. Rev. D_ 49, 1994, pp. 3209–3227  DOI: 10.1103/PhysRevD.49.3209
  * [27] Alexander L. Read  “Modified frequentist analysis of search results (The CL(s) method)”  In _Workshop on Confidence Limits_ , 2000, pp. 81–101
  * [28] Jerzy Neyman and Egon Sharpe Pearson  “On the Problem of the Most Efficient Tests of Statistical Hypotheses”  In _Phil. Trans. Roy. Soc. Lond. A_ 231.694-706, 1933, pp. 289–337  DOI: 10.1098/rsta.1933.0009
  * [29] LEP SUSY Working Group, ALEPH, DELPHI, L$3$ and OPAL Collaborations  “Stable Heavy Charged Particles”  URL: http://lepsusy.web.cern.ch/lepsusy/www/stable%5C_summer02/stable_208.html
  * [30] LEP SUSY Working Group, ALEPH, DELPHI, L$3$ and OPAL Collaborations  “Combined LEP GMSB Stau//Smuon/Selectron Results, 189-208 GeV”  URL: http://lepsusy.web.cern.ch/lepsusy/www/gmsb_summer02/lepgmsb.html
  * [31] Serguei Chatrchyan  “Searches for Long-Lived Charged Particles in $pp$ Collisions at $\sqrt{s}$=7 and 8 TeV”  In _JHEP_ 07, 2013, pp. 122  DOI: 10.1007/JHEP07(2013)122
  * [32] Andreas Haupt and Yves Kemp  “The NAF: National Analysis Facility at DESY”  In _J. Phys. Conf. Ser._ 219, 2010, pp. 052007  DOI: 10.1088/1742-6596/219/5/052007