# Deconvolution on graphs via linear programming

Bernhard G. Bodmann  641 Philip G. Hoffman Hall, Department of Mathematics,
University of Houston, Houston, TX 77204-3008<EMAIL_ADDRESS>and  Jennifer
J. May

###### Abstract.

The main challenge addressed in this paper is to identify individual terms in
a superposition of heat kernels on a graph. We establish geometric conditions
on the vertices at which these heat kernels are centered and find bounds on
the time parameter governing the evolution under the heat semigroup that
guarantee a successful recovery. This result can be viewed as a type of
deconvolution on a graph.

This research was supported in part by NSF Grant DMS ATD-1925352

##  1\. Introduction

Sparse recovery is generally understood as the identification of individual
terms in a linear combination of vectors from a dictionary. A main challenge
is to establish an accurate estimate of these terms even if one observes only
a noisy version of a low-dimensional projection of the given linear
combination. Under the theme of “compressed sensing off the grid” or
“superresolution”, this type of task has been accomplished for identifying
superpositions of sinusoids [17] or finitely supported measures from lowest
Fourier coefficients [4]. In yet another work, these ideas were developed in
the setting of Hilbert spaces with reproducing kernels [1] or in spline-type
spaces [18]. Here, we choose an application to graph signal processing [14,
15], where ideas from harmonic analysis are realized in a less structured
setting given by discrete geometries. In a few papers, the idea of sampling a
function and finding an accurate approximation is developed, mostly for a
class of smooth functions, see [13, 7, 6, 16]. Here, we focus on the recovery
of sparse functions, whose support is small, of which only a smoothed version
is observed. A main insight in this paper is that even in the discrete setting
of graphs, techniques from reproducing kernel spaces are useful.

In the main part of this paper, we assume that a function $g$ has support
$\\{v_{j}\\}_{j=1}^{J}$ and the values $g(v_{j})=c_{j}$,
$j\in\\{1,2,\dots,J\\}$. We wish to recover an approximation to $g$ from
observing the smoothed version $e^{t\Delta}g$ corrupted by noise,

$f=\sum_{j=1}^{J}c_{j}e^{t\Delta}\delta_{v_{j}}+w$

where $t\geq 0$ is the time parameter for the heat semigroup
$\\{e^{t\Delta}\\}_{t\geq 0}$ and $w$ is the additive noise. Applying a linear
inversion of $e^{t\Delta}$ is problematic. Even with a regularization step
using an orthogonal projection $P_{\Lambda}$, which projects on the sum of all
eigenspaces corresponding to eigenvalues of $-\Delta$ that are in the interval
$[0,\Lambda]$, the approximately recovered vector is

$\hat{g}=e^{-t\Delta}P_{\Lambda}f=\sum_{j=1}^{J}c_{j}P_{\Lambda}\delta_{v_{j}}+P_{\Lambda}e^{-t\Delta}w$

and the worst-case error norm is bounded below by
$|e^{t\Lambda}\|w\|-\|(I-P_{\Lambda})g\||$, which exhibits exponential growth
in $t$. A linear recovery procedure based on regularizing with $P_{\Lambda}$
and inverting the heat semigroup then needs to trade off errors between the
effect of $P_{\Lambda}$ on $g$ and the effect of the noise term $P_{\Lambda}w$
when the heat semigroup is inverted.

In this paper, we formulate a simple procedure that only requires an a priori
knowledge of the noise norm $\|w\|$ in order to find an accurate approximation
of the original sparse signal. To this end, we formulate a linear program that
minimizes the $\ell^{1}$-norm among all signals that are consistent with the
observed data. This strategy has an extensive history in applications of
compressed sensing [2, 3, 11].

######  Definition 1.

Let $V$ be a finite set, $f\in\ell^{2}(V)$ and $\epsilon\geq 0$. The vector
$\hat{g}\in\ell^{2}(V)$ is a solution of the linear program (LP) if
$\|\hat{g}\|_{1}\leq\|\tilde{g}\|_{1}$ for each $\tilde{g}$ with
$\|e^{t\Delta}\tilde{g}-f\|\leq\epsilon$.

This definition is chosen so that when $f=e^{t\Delta}g+w$ with
$g,w\in\ell^{2}(V)$ unknown and $\|w\|\leq\epsilon$, then $g$ is in the
feasible domain over which the $\ell^{1}$-norm is minimized.

We use dual certificates, elements in $\ell^{\infty}(V)$, to verify that a a
vector is a norm minimizer. To illustrate this, we first consider the
noiseless case.

######  Definition 2.

If $f,\tilde{g}\in\ell^{2}(V)$ satisfy $e^{t\Delta}\tilde{g}=f$ and there
exists $a\in\ell^{2}(V)$ such that $h=e^{t\Delta}a$ gives $\|h\|_{\infty}\leq
1$ and $\|\tilde{g}\|_{1}=\langle f,a\rangle$, then we say that $h$ is a dual
certificate for $\tilde{g}$.

In this case, $\tilde{g}$ is a solution of (LP) because
$\|\tilde{g}\|_{1}=\langle f,a\rangle=\langle g,h\rangle=|\langle
g,h\rangle|=\|g\|_{1}$ rules out that there is another feasible choice for
$\tilde{g}$ with a lower norm. Further below, we recall that the same argument
can also be used to show uniqueness and stability with respect to noise, when
the dual certificate satisfies additional conditions.

In our context, for $g=\sum_{j=1}^{J}c_{j}\delta_{v_{j}}$, we construct a dual
certificate in the form

$h(x)=\sum_{j=1}^{J}a_{j}e^{t\Delta}\delta_{v_{j}}\,.$

In the absence of noise, the conditions then imply $\|h\|_{\infty}=1$ and
$h(v_{j})=c_{j}/|c_{j}|$ for each $j\in\\{1,2,\dots,J\\}$. We argue
qualitatively why such a dual certificate exists for small time. If $t=0$,
then setting $h(x)=c_{j}/|c_{j}|$ if $x=v_{j}$ and otherwise $h(x)=0$ gives
the certificate. We denote $S=\\{v_{j}\\}_{j=1}^{J}$. The heat semigroup
$\\{e^{t\Delta}\\}_{t\geq 0}$ is invertible on $\ell^{2}(V)$, hence also on
$\ell^{2}(S)$. We denote this operator by $M^{t}$, and its inverse by
$M^{-t}$. To be precise, if $f\in\ell^{2}(S)$, then we extend $f$ by $f(x)=0$
if $x\not\in S$ and let for $y\in S$

$(M^{t}f)(y)=(e^{t\Delta}f)(y)\,.$

By continuity, $\|M^{-t}h\|_{\infty}$ is close to $\|h\|_{\infty}=1$ for all
sufficiently small $t$. Now extending $a\equiv M^{-t}h$ by zero on $V\setminus
S$ and using that $e^{t\Delta}a$ is continuous, for sufficiently small $t$,
$|(e^{t\Delta}a)(x)|<1$ for each $x\not\in S$. The remainder of this paper
makes this argument quantitative, with explicit conditions for the geometry of
the set $S$ and bounds on the time of the heat semigroup that permit the
construction of a dual certificate. The main resource for the results
presented here is controlling the decay of the heat kernel on the graph. With
bounds on the heat kernel and the discrete geometry of the graph, we have an
alternative to elements of harmonic analysis that have been used to construct
dual certificates [4].

The remainder of this paper is organized as follows: In Section 2, we fix
notation. In the interest of keeping the exposition accessible, we first treat
the noiseless case in Section 3, and then build the case for noisy recovery in
Section 4.

##  2\. Preliminaries

Much of the material under consideration is based on the concept of a weighted
graph and an associated graph Laplacian.

######  Definition 3.

An edge-weighted finite graph is a triple $(V,E,b)$ where $V$ is the finite
set of vertices, the edge set $E$ contains subsets of $V$ of size two, and
$b:V\times V\to[0,\infty)$ has the property $b(x,y)=b(y,x)>0$ if and only if
$\\{x,y\\}\in E$. We say that $d:V\times V\to[0,\infty)$ is a distance
function that is compatible with the weight $b$ if it is a metric on $V$ such
that for $x,y\in V$, $d(x,y)=\min_{z\in V}({d(x,z)+d(z,y)})$, and for each
$x\in V$,

$\sum_{y:\\{x,y\\}\in E}b(x,y)(d(x,y))^{2}\leq 1\,.$

This type of distance has been investigated by Davies in the study of
semigroups generated by graph Laplacians [9]. If, for example,
$b(x,y)=(d(x,y))^{-2}/k$ with $k$ the maximal degree of the graph, then the
definition of a graph Laplacian is given by a difference quotient with
similarities to Euclidean geometry.

######  Definition 4.

Let $(V,E,b)$ be an edge-weighted finite graph, then the associated graph
Laplacian $\Delta$ is defined in terms of the quadratic form
$f\mapsto\langle\Delta f,f\rangle$ on $\ell^{2}(V)$,

$\langle\Delta f,f\rangle=-\frac{1}{2}\sum_{\\{x,y\\}\in
E}b(x,y)|f(x)-f(y)|^{2}$

The operator $\Delta$ is Hermitian and hence we can define the associated heat
semigroup $\\{e^{t\Delta}\\}_{t\geq 0}$, either through the spectral
representation or by the power series expansion of the exponential function.
We also describe this semigroup in terms of its kernel function,

$K(t,x,y)\equiv(e^{t\Delta}\delta_{y})(x)\,.$

The problem we investigate is under which conditions we can identify $g$ when
observing $f=e^{t\Delta}g+w$.

This paper will show the quantitative argument for both the noiseless and
noisy cases.

##  3\. Recovery in the Noiseless Case

We claim for small $t$ and a minimum distance between support vertices that we
can create a dual certificate and thus recover the signal, $g$, from the
observed signal, $f$. The following result can be found in Candès and
Fernandez-Granda [4], as well as Tang [17].

######  Theorem 5.

Let $g\in\ell^{2}(V)$ with support $S=\\{x:g(x)\neq 0\\}$. Suppose there
exists $a\in\ell^{2}(V)$ such that $e^{t\Delta}a=:h$ has the following
properties:

  1. (1)

$h(v_{j})=c_{j}/|c_{j}|$ for all $v_{j}\in S$

  2. (2)

$|q(x)|<1$ for all $x\in V\backslash S$,

then $g$ is the unique solution to the $\ell^{1}$-minimization problem (LP)
satisfying $e^{t\Delta}g=f.$

###### Proof.

Assume such a $h$ exists and let $\tilde{g}$ be a solution to the minimization
problem described. Then, $\tilde{g}$ is such that
$e^{t\Delta}\tilde{g}=e^{t\Delta}g.$ Moreover,

(1) $\begin{split}\langle\tilde{g},h\rangle&=\langle
e^{t\Delta}\tilde{g},a\rangle\\\ &=\langle g,e^{t\Delta}a\rangle\\\ &=\langle
g,h\rangle.\\\ \end{split}$

Now, assuming $\tilde{g}\neq g$ then $\tilde{g}=g+w$, with $w\neq 0$ and
$\langle w,h\rangle=0$. We can now decompose $w$ into a part that is zero
outside of the support of $g$ and the remainder, $w=w_{S}+w_{S^{c}}$. Then

$0=\langle h,w\rangle=\langle h,w_{S}\rangle+\langle h,w_{S^{c}}\rangle\,,$

so

$|\langle h,w_{S}\rangle|=|\langle h,w_{S^{c}}\rangle|$

and hence

$\lVert w_{S}\rVert_{1}=|\langle h,w_{S^{c}}\rangle|\leq\lVert
w_{S^{c}}\rVert_{1}\cdot\lVert h\big{|}_{S^{c}}\rVert_{\infty}.$

Using property (2) for $h$ then

$\lVert w_{S}\rVert_{1}<\lVert w_{S^{c}}\rVert_{1}.$

Now considering

(2) $\begin{split}\lVert g\rVert_{1}&\geq\lVert\tilde{g}\rVert_{1}=\lVert
g+w\rVert_{1}\\\ &=\lVert g+w_{S}+w_{S^{c}}\rVert_{1}\\\ &=\lVert
g+w_{S}\rVert_{1}+\lVert w_{S^{c}}\rVert_{1}\\\ &\geq\lVert g\rVert_{1}-\lVert
w_{S}\rVert_{1}+\lVert w_{S^{c}}\rVert_{1}\\\ &>\lVert g\rVert_{1}\\\
\end{split}$

gives a contradiction and hence $w=0$, proving the solution is unique. ∎

The next step is to establish the existence of such a dual certificate. We
recall that we assume a distance function on the graph that is compatible with
the weight in the definition of the graph Laplacian.

To prepare the construction of the dual certificate, we derive a lower bound
for the diagonal entries of the heat kernel based on spectral properties of
the Laplacian.

######  Lemma 6.

Suppose $G=(V,E,b)$ is an edge-weighted finite connected graph that has $N$
vertices and a subset $S\subset V$ with $S=\\{v_{1},\dots,v_{J}\\}$ for some
$J\in\mathbb{N}$. Let $\Delta$ be the graph Laplacian of $G$ with operator
norm $\|\Delta\|$ and let $\lambda$ be the smallest non-zero eigenvalue of
$-\Delta$, then for $x\in V$,

$\frac{1}{N}+e^{-t\lVert\Delta\rVert}\left(\frac{N-1}{N}\right)\leq
K\left(t,x,x\right)\leq\frac{1}{N}+e^{-t\lambda}\left(\frac{N-1}{N}\right)\,.$

###### Proof.

We recall that the constant function is in the eigenspace of $\Delta$
corresponding to eigenvalue zero. We can split $\delta_{x}$ in a linear
combination of a constant and a zero-summing vector.

  1. (1)

Let $\mu(x)=1/N$, $\nu_{x}=\delta_{x}-\mu$, then $\delta_{x}=\mu+\nu_{x}$, and
$\Delta\mu=0$.

  2. (2)

Using the orthogonality of $\mu$ and $\nu_{x}$, we compute the (squared)
norms, yielding $\|\mu\|^{2}=\frac{1}{N}$ and

$\lVert\nu_{x}\rVert^{2}=1-\frac{1}{N}=\frac{N-1}{N}\,.$

Next, we consider the diagonal entries of the heat semigroup:

(3) $\displaystyle K(t,x,x)$ $\displaystyle=\langle
e^{t\Delta}\delta_{x},\delta_{x}\rangle$ (4) $\displaystyle=\langle
e^{t\Delta}(\mu+\nu_{x}),(\mu+\nu_{x})\rangle$ (5) $\displaystyle=\langle
e^{t\Delta}\mu,\mu\rangle+\langle e^{t\Delta}\mu,\nu_{x}\rangle+\langle
e^{t\Delta}\nu_{x},\mu\rangle+\langle e^{t\Delta}\nu_{x},\nu_{x}\rangle$ (6)
$\displaystyle=\langle\mu,\mu\rangle+\langle
e^{t\Delta}\nu_{x},\nu_{x}\rangle$ (7)
$\displaystyle\geq\lVert\mu\rVert^{2}+e^{-t\lVert\Delta\rVert}\lVert\nu_{x}\rVert^{2}$
(8)
$\displaystyle=\frac{1}{N}+e^{-t\lVert\Delta\rVert}\left(\frac{N-1}{N}\right)\,.$

A similar line of reasoning gives

$K(t,x,x)\leq\frac{1}{N}+e^{-t\lambda}\left(\frac{N-1}{N}\right)\,.$

∎

We rely on a work by Folz [10], building on Davies [8] to bound the off
diagonal entries of the heat kernel matrix.

######  Lemma 7 (Theorem 2.1 [10]).

Let $(V,E,b)$ be an edge-weighted finite graph and $d$ a compatible distance
function, then the heat kernel has the bound for each $x,y\in V$ and $t>0$,

$K(t,x,y)\leq\exp\Big{[}-\frac{1}{2}d(x,y)\log\left(\frac{d(x,y)}{2et}\right)\Big{]}\,.$

The bounds on the heat kernel establish invertibility of the operator on
$\ell^{2}(S)$, if the vertices in $S$ are sufficiently separated.

######  Proposition 8.

Given an edge-weighted finite connected graph $G=(V,E,b)$ with a compatible
distance function $d$, the graph Laplacian $\Delta$ and the operator $M^{t}$
associated with the heat kernel on $\ell^{2}(S)$, if $|S|=J$, and for each
$\\{x,y\\}\subset S$ we have $d(x,y)\geq D$, and if $t\leq T$ with

$(J-1)\left(\frac{2eT}{D}\right)^{D/2}<\frac{1}{N}+e^{-T\lVert\Delta\rVert}\left(\frac{N-1}{N}\right),$

then $M^{t}$ is invertible on $\ell^{2}(S)$ and the operator norm of the
inverse is bounded by

$\|M^{-t}\|\leq\Bigl{(}\frac{1}{N}+e^{-T\lVert\Delta\rVert}\left(\frac{N-1}{N}\right)-(J-1)\left(\frac{2eT}{D}\right)^{D/2}\Bigr{)}^{-1}\,.$

###### Proof.

We note that the inequality implies that $D>2eT$, otherwise this contradicts
with the right-hand side being bounded above by one. Thus, we can use
monotonicity in $d(x,y)\geq D$ to further estimate Folz’s result,

$K(t,x,y)\leq\exp\left[-\frac{1}{2}D\log\left(\frac{D}{2et}\right)\right]=\left(\frac{2et}{D}\right)^{D/2}$

and then use monotonicity in $t$ to bound off-diagonal elements of $M^{t}$ for
each $t\leq T$. Now using the Gersgorin’s theorem [12, Thm 6.2.27], $M^{t}$ is
invertible because it is symmetric and each eigenvalue is bounded below by the
minimum among the diagonal entries minus the largest sum occurring among the
$(J-1)$ off-diagonal entries in each row. Consequently, the operator norm of
the inverse is bounded by the inverse of the lower bound for the smallest
eigenvalue of $M^{t}$. ∎

Now that we have shown the heat semigroup is invertible on the vertices in the
support of our signal, we look towards including the non-support vertices. To
do so, we will use a classical result from Varah [19], which controls the
$\ell^{\infty}$ norm under $M^{-t}$.

######  Lemma 9.

Given $G$, $\Delta$, $S$ and $M^{t}$ as above, then if
$h(v_{j})=c_{j}/|c_{j}|$ for $j\in\\{1,2,\dots,J\\}$ and $a=M^{-t}h$, we have

$\|a\|_{\infty}\leq\Bigl{(}\frac{1}{N}+e^{-T\|\Delta\|}\frac{N-1}{N}-(J-1)(\frac{2eT}{D})^{D/2}\Bigr{)}^{-1}$

###### Proof.

If a matrix $A=(a_{k,j})_{k,j=1}^{J}$ is strictly diagonally dominant matrix,
then the operator norm of its inverse on $\ell^{\infty}$ is bounded by

$\lVert A^{-1}\rVert_{\infty,\infty}<\left(\min_{k}\left(\lvert
a_{kk}\rvert-\sum_{j\neq k}\lvert a_{kj}\rvert\right)\right)^{-1}.$

Now inserting the heat semigroup $M^{t}$ for $A$ gives in our case that if
$|h||_{S}=1$, then $a=M^{-t}h$ has

$\|a\|_{\infty}\leq\bigl{(}\min_{x}(K(t,x,x)-(J-1)\sum_{y:y\neq
x}K(t,x,y)\bigr{)}^{-1}\,.$

With Folz’s bound, this yields again the same explicit expression as for the
$\ell^{2}$-norm,

$\|a\|_{\infty}\leq\Bigl{(}\frac{1}{N}+e^{-T\|\Delta\|}\frac{N-1}{N}-(J-1)\bigl{(}\frac{2eT}{D}\bigr{)}^{D/2}\Bigr{)}^{-1}\,.$

∎

######  Theorem 10.

Suppose we have a finite edge-weighted connected graph, $G=(V,E,b)$ and a
compatible distance function. Let $\Delta$ be the graph Laplacian of $G$ and
$N=|V|$ Given a subset of vertices of size $J$ with a minimum distance $D$
between any two vertices in

$S:=\\{v_{1},\dots,v_{J}\\}\text{with }d(v_{j},v_{k})\geq D,\mathrm{if}j\neq
k$

and let $\zeta$ be the smallest distance between vertices in $V$.

If the following relations hold for $T$ and $D$,

  1. (1)

$(J-1)\left(\frac{2eT}{D}\right)^{D/2}<\frac{1}{N}+e^{-T\lVert\Delta\rVert}\left(\frac{N-1}{N}\right)$
and

  2. (2)

$\left(\frac{2eT}{\zeta}\right)^{\zeta/2}+\left(\frac{4eT}{D}\right)^{D/4}\cdot(J-1)<\frac{1}{N}+e^{-T\lVert\Delta\rVert}\left(\frac{N-1}{N}\right)-(J-1)\left(\frac{2eT}{D}\right)^{D/2}$,

then for any choice $\epsilon_{j}\in\\{\pm 1\\}$, $j\in\\{1,2,\dots,J\\}$,
there is a linear combination,
$h(x)=\displaystyle\sum_{j=1}^{J}a_{j}e^{t\Delta}(x,v_{j})$, with the
following properties:

  1. (i)

$\lVert h\rVert_{\infty}=1$

  2. (ii)

$h(v_{j})=\epsilon_{j}$, for all $v_{j}\in S$

  3. (iii)

$\lvert h(x)\rvert<1$ for all $x\in V\backslash S$.

###### Proof.

Note: As defined above, $h_{j}=\left(M^{t}a\right)_{j}$ and $\lVert
h\rVert_{\infty}=1.$ Using Proposition 8, $M^{t}$ is strictly diagonally
dominant and Lemma 9 applies, with $a=M^{-t}h$. Thus

$\lVert
a\rVert_{\infty}=\lVert{M}^{-t}h\rVert_{\infty}\leq\lVert{M}^{-t}\rVert_{\infty}\cdot\lVert
h\rVert_{\infty}=\lVert{M}^{-t}\rVert_{\infty}.$

Next, we extend the heat kernel from $M^{t}$ to $K\left(t,x,y\right)$, where
$x,y\in V$. We wish to show $\lvert h(v)\rvert<1$, for every $v\notin S$,
knowing we have picked $a_{j}$’s that make $h(v_{j})\in\\{-1,1\\}$ for all
$v_{j}\in S.$

So consider an arbitrary $v\not\in S.$ By the separation of elements in $S$,
there exists at most one $v_{1}\in S$ such that $\zeta\leq d(v,v_{1})\leq D/2$
and for all other $v_{j}\in S$ we have $d(v,v_{j})\geq D/2.$ Using Folz’s
bound [10] we have

(10)
$K(t,v,v_{1})\leq\left(\frac{2et}{\zeta}\right)^{\zeta/2}\text{\leavevmode\nobreak\
and \leavevmode\nobreak\ \leavevmode\nobreak\
}K(t,v,v_{j})\leq\left(\frac{4et}{D}\right)^{D/4},$

and if $v\not\in S$,

(11) $\displaystyle\lvert h(v)\rvert$ $\displaystyle<\lVert
a\rVert_{\infty}\left(\left(\frac{2et}{\zeta}\right)^{\zeta/2}+(J-1)\left(\frac{4et}{D}\right)^{D/4}\right)<1.$

∎

######  Corollary 11.

Under the assumptions of the preceding theorem, a solution $g$ to (LP) is the
unique function $g$ such that $f=e^{t\Delta}g$.

##  4\. Recovery with Noisy Measurements

Next, we consider the effect of noise. It is fortuitous that there is no need
to alter the construction of the dual certificate in the treatment of a noisy
measurement. We proceed as described in the Introduction and find a solution
to (LP). If the noise norm is known, so

$f=e^{t\Delta}g+w$

has been observed with $g,w$ unknown but $\|w\|\leq\epsilon$ then the accuracy
of the recovered vector is controlled proportional to $\epsilon$. We will show
that one can stably recover $g$ by finding a $\hat{g}$ minimizer to (LP),
which then satisfies $\|\hat{g}-g\|\leq C\epsilon$ with some fixed constant
$C$.

We begin with a lemma. To simplify notation, we abbreviate
$\|M^{-t}\|=1+\delta$. We also split
$\ell^{2}(V)=\ell^{2}(S)\oplus\ell^{2}(S^{c})$ and denote the corresponding
terms $f=f_{S}+f_{S^{c}}$, where the index contains the support of each term.

######  Lemma 12.

Let $f=e^{t\Delta}g+w$ with $\|w\|\leq\epsilon$ and $\hat{g}$ be a minimizer
of (LP), then

$\|\hat{g}-g\|_{1}\leq
2\sqrt{J}(1+\delta)\epsilon+\|{\hat{g}}_{S^{c}}\|_{1}\,.$

###### Proof.

By the splitting and the fact that $S$ contains the support of $g$,

$\|\hat{g}-g\|_{1}=\|\hat{g}_{S}-g_{S}\|_{1}+\|{\hat{g}}_{S^{c}}\|_{1}\,.$

We further estimate the first term using the Cauchy-Schwarz inequality,

$\|\hat{g}_{S}-g_{S}\|_{1}\leq\sqrt{J}\|\hat{g}_{S}-g_{S}\|_{2}\leq\sqrt{J}\|M^{-t}\|\|M^{t}(\hat{g}_{S}-g_{S})\|\leq\sqrt{J}(1+\delta)\|e^{t\Delta}(\hat{g}-g)\|\leq
2\sqrt{J}(1+\delta)\epsilon\,.$

The last inequality is due to the tolerance in the feasible region and the
second last inequality follows from $M^{t}=\chi_{S}e^{t\Delta}\chi_{S}$, where
$\chi_{S}$ is the multiplication operator with the characteristic function of
$S$, projecting onto $\ell^{2}(S)$. The operator inequality
$\chi_{S}e^{t\Delta}\chi_{S}\leq e^{t\Delta}$ then implies
$\|M^{t}\eta_{S}\|\leq\|e^{t\Delta}\eta\|$ for any $\eta\in\ell^{2}(V)$. ∎

Henceforth, we abbreviate $\eta=g-\hat{g}$. The following lemmata are similar
to results in Candès and Plan [5]. However, we keep using $h$ and do not
construct an approximate dual certificate, even in the noisy case.

######  Lemma 13.

With a dual certificate $h$ for $g$, the following norm inequality holds:

(12) $\lVert\hat{g}\rVert_{1}\geq\lVert
g\rVert_{1}-\langle\eta_{S},h_{S}\rangle+\|\hat{g}_{S^{c}}\|_{1}.$

###### Proof.

We define $\hat{h}(x)=\hat{g}(x)/|\hat{g}(x)|$ for
$x\in\mathrm{supp}(\hat{g})$ and otherwise $\hat{h}(x)=0$, then

(13)
$\begin{split}\lVert\hat{g}\rVert_{1}&=\langle\hat{g}_{S},\hat{h}_{S}\rangle+\|\hat{g}_{{S^{c}}}\|_{1}\\\
&\geq\langle\hat{g}_{S},h_{S}\rangle+\lVert\hat{g}_{S^{c}}\rVert_{1}+\langle
g_{S},h_{S}\rangle-\langle g_{S},h_{S}\rangle\\\ &=\lVert
g\rVert_{1}-\langle\eta_{S},h_{S}\rangle+\lVert\hat{g}_{S^{c}}\rVert_{1}.\end{split}$

∎

######  Lemma 14.

If $\hat{g}$ is a solution of (LP), then

$\|\hat{g}_{S^{c}}\|_{1}\leq|\langle\eta_{S},h_{S}\rangle|\,.$

###### Proof.

This is because $\hat{g}$ minimizes the $\ell^{1}$-norm in the feasible set,
so $\|\hat{g}\|_{1}\leq\|g\|_{1}$. Consequently, the previous lemma implies

$0\geq\|\hat{g}\|_{1}-\|g\|_{1}\geq-\langle\eta_{S},h_{S}\rangle+\lVert\hat{g}_{S^{c}}\rVert_{1},$

which simplifies to the claimed inequality. ∎

We are ready to conclude the error bound for sparse recovery from noisy
measurements.

######  Theorem 15.

Given a finite edge-weighted connected graph $G=(V,E,b)$, and a $J$-sparse
signal $g$ on $G$. Assume $T$, $J$, and $D$ are as in Theorem 10 and $\hat{g}$
is a solution to the linear program (LP), then

(14) $\|\hat{g}-g\|_{2}\leq\|\hat{g}-g\|_{1}\leq
4(1+\delta)\sqrt{J}\epsilon\,.$

###### Proof.

Combining the lemmas, we already know

$\|\hat{g}-g\|_{1}\leq
2\sqrt{J}(1+\delta)\epsilon+|\langle\eta_{S},h_{S}\rangle|\,.$

We consider the last term,

$|\langle\eta_{S},h_{S}\rangle|=|\langle M^{t}\eta_{S},a\rangle|$

where $h|_{S}=M^{t}a$. Now using the Cauchy-Schwarz inequality and again
extending $M^{t}\eta_{S}$ to $e^{t\Delta}\eta$ on all of $V$ gives

$|\langle\eta_{S},h_{S}\rangle|\leq\|\sqrt{J}e^{t\Delta}\eta\|\|a\|_{\infty}\leq
2\epsilon\sqrt{J}(1+\delta)\,.$

∎

## References

  * [1] B. G. Bodmann, A. Flinth, and G. Kutyniok.  Compressed sensing for analog signals, https://arxiv.org/abs/1803.04218, 2018.
  * [2] E. Candès, J. Romberg, and T. Tao.  Stable signal recovery from incomplete and inaccurate measurements.  Commun. Pure Appl. Math., pages 1207–1223, 2005.
  * [3] E. Candès, M. Rudelson, T. Tao, and R. Vershynin.  Error correction via linear programming.  In 46th Annual IEEE Symposium on Foundations of Computer Science (FOCS’05), pages 668–681. IEEE, 2005.
  * [4] E. J. Candès and C. Fernandez-Granda.  Towards a mathematical theory of super-resolution.  Communications on Pure and Applied Mathematics, 67(6):906–956, 2014\.
  * [5] E. J. Candès and Y. Plan.  A probabilistic and ripless theory of compressed sensing.  IEEE transactions on Information Theory, 57(11):7235–7254, 2011\.
  * [6] L. F. O. Chamon and A. Ribeiro.  Greedy sampling of graph signals.  IEEE Trans. Signal Process., 66(1):34–47, 2018.
  * [7] S. Chen, A. Sandryhaila, J. M. F. Moura, and J. Kovačević.  Signal recovery on graphs: variation minimization.  IEEE Trans. Signal Process., 63(17):4609–4624, 2015.
  * [8] E. B. Davies.  Large Deviations for Heat Kernels on Graphs.  Journal of the London Mathematical Society, s2-47(1):65–72, 02 1993\.
  * [9] E.B. Davies.  Analysis on graphs and noncommutative geometry.  Journal of functional analysis, 111(2):398–430, 1993.
  * [10] M. Folz.  Gaussian upper bounds for heat kernels of continuous time simple random walks.  Elec. J. Prob., 62:1693–1722, 2011.
  * [11] S. Foucart and H. Rauhut.  A Mathematical Introduction to Compressive Sensing.  Birkhauser, Boston, 2013.
  * [12] R. A. Horn and C. R. Johnson.  Matrix analysis.  Cambridge University Press, Cambridge [Cambridgeshire] ;, 1985.
  * [13] I. Z. Pesenson and M. Z. Pesenson.  Sampling, filtering and sparse approximations on combinatorial graphs.  J. Fourier Anal. Appl., 16(6):921–942, 2010.
  * [14] A. Sandryhaila and J. M. F. Moura.  Discrete signal processing on graphs.  IEEE Trans. Signal Process., 61(7):1644–1656, 2013.
  * [15] A. Sandryhaila and J. M. F. Moura.  Discrete signal processing on graphs: frequency analysis.  IEEE Trans. Signal Process., 62(12):3042–3054, 2014.
  * [16] Y. Tanaka.  Spectral domain sampling of graph signals.  IEEE Trans. Signal Process., 66(14):3752–3767, 2018.
  * [17] G. Tang, B. N. Bhaskar, P. Shah, and B. Recht.  Compressed sensing off the grid.  IEEE transactions on information theory, 59(11):7465–7490, 2013\.
  * [18] M. Unser.  A unifying representer theorem for inverse problems and machine learning.  Foundations of Computational Mathematics, 21(4):941–960, 2021.
  * [19] J.M. Varah.  A lower bound for the smallest singular value of a matrix.  Linear Algebra and its Applications, 11(1):3–5, 1975.


# Regression Test Suite for Payment Switch using jPOS

A. Sardesai, H. Piplodwala, S. Hargunani, S. Sonawane, Prof. P. Joshi, V.
Mohite  Computer Engineering Department, Pune Institute of Computer Technology

###### Abstract

The Payment Switch is an integral component of all modern payment and banking
systems in India. The NPCI currently provides a simulator to test payment
switches. However, this system has a few disadvantages viz. it lacks an API,
it requires manual generation of each test case and during high server loads,
the testing process may take a long time. Currently there aren’t any open
source alternatives to the NPCI simulator. We propose a system which solves
these shortcomings. Our proposed system simulates the NPCI system. It allows
connection with switches that are to be tested and automates the process of
generation and execution of test cases. It also has the capability to generate
test reports and can be run locally.

###### Index Terms:

Financial Technology, Regression Testing, Quality Assurance, Automation

##  I Introduction

Financial technology is rapidly advancing and with it, payments have become
synonymous with credit and debit card transactions. These transactions are
being fuelled by modern payment systems. The switch is an integral and
inseparable component of these systems. The NPCI (National Payments
Corporation of India) currently uses the ISO8583 platform to relay messages
between payment switches. The Switch generally offers a rules-based
authorisation and switching solution that is merchant driven. Payment
transactions are constantly routed between numerous acquirers and Payment
Service Providers. The switch is critical in the processing and validation of
various types of payments.

Such an important component would naturally require frequent updates from
switch developers. Every update must be tested rigorously; otherwise it might
result in invalid transactions and immense losses for the switch provider.
Currently, evaluating these switches is a time-consuming and laborious
process. It entails hundreds of test cases being conducted by individual
switch developers. The NPCI has provided a simulator which allows testing of
bank switches but it has the following problems:

  1. 1.

It offers a GUI which requires manual input of test cases from the tester.

  2. 2.

It does not have an API which can allow automation of test case execution.

  3. 3.

During high server loads or bad internet connectivity, the tests can be very
slow to execute, increasing the overall testing time.

  4. 4.

When the simulator is offline, testing cannot be performed.

The ISO8583 protocol and building payment switches requires niche knowledge.
Hence, an open source implementation to solve this problem currently does not
exist.

The jPOS library in Java is an open source library which allows us to create
applications that deal with the ISO8583 protocol. We propose building a jPOS
Q2 application that solves the aforementioned problems.

##  II Related Work

In terms of work and computational resources, in [1] Prabhakar K. et al.
discuss how re-testing is costly for maintenance. As a result, this study
provides a comprehensive strategy to generating cases for regression testing,
which predicts test work needed and noticing faults.

In [2] Prof. A. Ananda Rao and Kiran Kumar J delve deeper into cost-effective
regression testing. The goal of this research is to figure out how to cut the
cost of regression testing. In a black box environment, a method for reducing
test suites for regression testing has been proposed.

In [3] Susanne Rösch, Sebastian Ulewicz, Julien Provost and Birgit Vogel-
Heuser examined and contrasted the various testing methodologies utilised in
production automation. The purpose is to evaluate these approaches in terms of
their applicability to the field of industrial automation in order to spot
current trends and research needs.

In [4] B. Korel, L.H. Tahat and B. Vaysburg proposed a model-based regression
testing strategy that reduces regression test suites by using EFSM (Extended
Finite State Machine) model dependence analysis.

In [5] Yanping Chen, Robert L. Probert and Hasan Ural discuss extended
dependence analysis. Based on EFSM dependence analysis, a model-based
regression test suite (RTS) reduction strategy is suggested.

##  III Overview

We propose building a jPOS Q2 application that solves the aforementioned
problems. The general attributes of this application can be stated as follows:

  1. 1.

It simulates the NPCI switch. (national financial switch not NPCI)

  2. 2.

It is able to connect to other switches.

  3. 3.

It can send and receive ISO-Messages using various channels like ASCII, NAC
and XML.

  4. 4.

Allows testers to define test case templates and automates test case
generation.

  5. 5.

Generates a test report based on the execution of test cases after comparing
the received response against the expected response.

  6. 6.

Allows users to interact with the system via a GUI or CLI.

  7. 7.

Can operate as a Server or Client based on requirements.

  8. 8.

Can be run locally and an external internet connection is not required.

Figure 1: High Level Design

Further details are provided in the remainder of this paper for specific
modules.

##  IV System Modules

We propose a system with the following key modules:

###  IV-A Message Generator

Testers need to create a template for a particular test case. This is defined
in a simple json file. The test case denotes an ISO8583 message. The randomize
attribute allows us to randomize certain fields.

The module is responsible for generating a message according to the ISO8583
specification. The message generator will use the json template to generate an
ISOMsg instance while randomizing the necessary fields. A default field config
file decides the regex for each field that needs to be randomized. This can be
customized according to your custom packager.

Figure 2: Message Generator - template for a sample transaction

###  IV-B Server

The Server is a jpos Q2 application. It consists of three servers running on
different ports. Each server caters to a specific ISO Channel (XML, ASCII,
NAC). The transaction manager has a custom participant that validates the
requests and returns an appropriate response. This is also where all of the
business logic for request handling is implemented.

A function inside our custom participant to validate balance inquiry
transactions:

procedure validateBalanceEnquiry()
1 let response := ”00”

2 if request[0] != ”0200”:
2.1 response := ”12”

3 if request[3] = ”31[0-9]4”:
3.1 response := ”12”

4 request.setMTI(”0210”)

5 request[39] := response

6 if response = ”00”
6.1 request[54] := balance

7 return request

###  IV-C Client

This is a jPOS Q2 application that is responsible for carrying out regression
testing on the payment switch. The Client allows the user to select multiple
channels to send transactions, as well as the number of message iterations. It
asynchronously sends all the test transactions to the switch simultaneously
and generates a report based on the switch’s responses. The key functions
performed by the client are:

  1. 1.

Establish connection with switch

  2. 2.

Pack test transaction in ISO format

  3. 3.

Decode response and compare with expected output

  4. 4.

Generate report based on responses seen

Figure 3: Client sending a transaction request to the server

Fig. 3 shows the client sending a transaction request to the server, indicated
by the isomsg direction value ”outgoing”. The request is sent through the
ASCII channel on port 8001. The request body contains important information
pertaining to the specific type of request, like terminal ID, MTI, bank
branch, card number etc.

###  IV-D User Interface System

We added a simple user interface system using Java Swing to ensure ease of
usage for switch developers. This allows switch developers to perform
regression testing on their system with just a few clicks. It also displays
reports and allows them to be downloaded for future reference. The main
features can be classified as:

  1. 1.

Allows selecting whether simulator should act as client or server

  2. 2.

Call client and server modules

  3. 3.

Display and allow downloading of reports

Figure 4: Client GUI Figure 5: Report Generation GUI

##  V Conclusion

Our proposed system can be used to extensively test both key roles of a
payment switch, viz. when the bank is an issuer bank as well as when the bank
is the acquirer. This was our key objective and we’re pleased to have met it.

Given that we send messages asynchronously, we have massively sped up the
process of testing. This improves switch developer productivity as they no
longer need to spend hours upon hours manually testing the switch for
potential bugs. All in all, we’re excited to present a system we built that
can test such an integral component of payment systems.

## References

  * [1] S. Kadry and A. Kalakech, “Cost-effectiveness of regression testing: Problem and solution,” vol. 6, pp. 328–335, 05 2011.
  * [2] A. Rao and K. Jogu, “An approach to cost effective regression testing in black-box testing environment,” _International Journal of Computer Science Issues_ , vol. 8, 05 2011.
  * [3] S. Rösch, S. Ulewicz, J. Provost, and B. Vogel-Heuser, “Review of model-based testing approaches in production automation and adjacent domains—current challenges and research gaps,” _Journal of Software Engineering and Applications_ , vol. 08, pp. 499–519, 01 2015.
  * [4] B. Korel, L. Tahat, and B. Vaysburg, “Model based regression test reduction using dependence analysis,” 10 2002, pp. 214–.
  * [5] Y. Chen, R. L. Probert, and H. Ural, “Regression test suite reduction using extended dependence analysis,” in _Fourth International Workshop on Software Quality Assurance: In Conjunction with the 6th ESEC/FSE Joint Meeting_ , ser. SOQUA ’07. New York, NY, USA: Association for Computing Machinery, 2007, p. 62–69. [Online]. Available: https://doi.org/10.1145/1295074.1295086


# Directional direct detection of light dark matter up-scattered by cosmic
rays from direction of the Galactic center

Keiko I. Nagao  Satoshi Higashino  Tatsuhiro Naka  Kentaro Miuchi

###### Abstract

Dark matter with MeV scale mass is difficult to detect with standard direct
search detectors. However, they can be searched for by considering the up-
scattering of kinetic energies by cosmic rays. Because the dark matter density
is higher in the central region of the Galaxy, the up-scattered dark matter
will arrive at Earth from the direction of the Galactic center. Once the dark
matter is detected, we can expect to recognize this feature by directional
direct detection experiments. In this study, we simulate the nuclear recoils
of the up-scattered dark matter and quantitatively reveal that a large amount
of this type of dark matter is arriving from the direction of the Galactic
center. Also, we have shown that the characteristic signatures of the up-
scattered dark matter can be verified with more than 5$\sigma$ confidence
levels for the assumed target atoms and future upgrades to directional
detectors.

##  1 Introduction

The identity of dark matter (DM) is a highly significant question in cosmology
and particle physics. In the past a few decades, weakly interacting massive
particles (WIMPs) have been expected to be the identity of DM. Many attempts
have been made to detect WIMPs directly and indirectly, and to produce them in
collider experiments. However, despite advances in the sensitivity of direct
detection experiments, no obvious detection of DM has been achieved so far.
Particularly, direct detection experiments have placed very stringent limits
on the magnitude of the interaction cross section for WIMP particles of GeV-
scale mass, i.e., $\sigma_{\mathrm{SI}}<6.5\times 10^{-48}$ cm2 and
$\sigma_{\mathrm{SD}}<3.1\times 10^{-41}$ cm2, where
$\sigma_{\mathrm{SI(SD)}}$ represents spin-independent (SI) and spin-dependent
(SD) cross sections of the DM and proton scattering, respectively [1, 2]. On
the other hand, for DM having mass lighter than the GeV-scale, it is not easy
to obtain high sensitivity through DM-nucleon scatterings. It has recently
been proposed that even in such light DM cases, unavoidable scattering between
DM and a cosmic ray produces substantial and very energetic DM that can be
detected in neutrino detectors and direct detection experiments of DM [3, 4,
5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. As
well as the up-scattering of DM by cosmic rays, its effect on cosmic ray
propagation in the Earth [6] has also been investigated. The cosmic ray
scattered dark matter (CR-DM) flux has a substantial contribution from the
direction of the Galactic center [6, 9, 18]. In some direct searches such as
PandaX and CDEX, constraints for CR-DM have already been placed using real
data [16, 17].

Directional direct detection of DM is expected to be very well suited to
detect the characteristic events of CR-DM coming from the Galactic center.
Using directional information to understand the DM nature has been
investigated in both directional and non-directional direct searches [24, 25,
26, 27, 28, 29, 30]. Also, verifications of CR-DM using directional
information through diurnal modulation in the ordinary direct detection
experiments [6, 16], and directional search in the neutrino detector [19] are
already studied. In this study, we simulate the CR-DM signals in directional
searches and reveals its potential for validation of CR-DM. In directional
detection, gas detectors are sensitive to the SD interactions between DM
particles and nucleons [31, 32]. The gases used in the experiments are CF4 and
SF6, and in most cases, the nucleon target for scattering with DM is fluorine
(F). We also simulate a solid detector NEWSdm [33, 34], which is a nuclear
emulsion detector and is sensitive to SI interactions. Nuclear emulsion
detectors contain multiple targets, thus in the numerical calculations we
employ hydrogen ($p$), which is sensitive to light DM, and silver (Ag), which
is sensitive to heavy DM.

The rest of this paper is organized as follows. Section 2 introduces the
methods of the CR-DM simulations and three assumed DM profiles. Section 3
presents numerical results and discusses the detectability. Finally, we
conclude the paper in Section 4.

##  2 CR-DM and the density profile

To simulate the scattering of CR-DM and the targets in direct detection
experiments, first a parametrization of the incoming flux of CR-DM to the
Earth is required. In this section, we follow the methods of reference [4] to
calculate the incoming flux of CR-DM from all directions in space. The
calculation procedure is listed below.

  1. 1.

If we assume the DM mass $m_{\chi}$, the interaction cross section
$\sigma_{\chi i}$ between DM and the nucleus $i$, density profile
$\rho_{\chi}(r)$, and cosmic ray flux $d\Phi_{i}/dT_{i}$ in the Galaxy, the
incoming CR-DM flux can be obtained by

$\frac{d\Phi_{\chi}}{dT_{\chi}d\Omega}=\int_{l.o.s.}l^{2}dl\frac{1}{4\pi
l^{2}}\sigma_{\chi
i}\frac{\rho_{\chi}(r)}{m_{\chi}}\int_{T_{i}^{\mathrm{min}}}^{\infty}dT_{i}\frac{d{\Phi_{i}}}{dT_{i}},$
(2.1)

where $l$ is the length of the line-of-sight (l.o.s) integral, $r$ is the
distance from the Galactic center. $T_{\chi}$ and $T_{i}$ are the kinetic
energies of the CR-DM and nucleus, respectively. Minimum kinetic energy
$T_{i}^{\mathrm{min}}$ is defined as

$T_{i}^{\mathrm{min}}=\left(\frac{T_{\chi}}{2}-m_{i}\right)\left(1\pm\sqrt{1+\frac{2T_{\chi}}{m_{\chi}}\frac{(m_{i}+m_{\chi})^{2}}{(2m_{i}-T_{\chi})^{2}}}\right)$

where $+$ ($-$) corresponds to the case for $T_{\chi}>2m_{i}$
($T_{\chi}<2m_{i}$). Additionally, we use Hydrogen, which is the most abundant
element in space, to be the cosmic ray particle $i$, i.e., we assume $i=p$ for
simplicity in this paper.

To define the direction of the incoming CR-DM and the nuclear recoils, we
adopt the coordinate system shown in Figure 1. We define the direction of the
incoming CR-DM with $\theta_{\mathrm{DM}}$ and $\phi_{\mathrm{DM}}$, and the
direction of the nuclear recoils as $\theta_{\mathrm{N}}$ and
$\phi_{\mathrm{N}}$. In this coordinate system,
$\theta_{\mathrm{DM}}=\phi_{\mathrm{DM}}=0$ corresponds to events coming from
the direction of the Galactic center. The scattering angle
$\gamma_{\mathrm{G.C.}}$ is the angle from the line connecting the Galactic
center to the observer. There is a relation
$\cos{\gamma_{\mathrm{G.C.}}=\cos{\theta_{\mathrm{N}}}\cos{\phi_{\mathrm{N}}}}$
between $\theta_{\mathrm{N}},\phi_{\mathrm{N}}$ and $\gamma_{\mathrm{G.C.}}$.
To include contributions near the Galactic center in the integration of $l$,
we take the limits of integration up to the distance between the solar system
and the Galactic center $d$ and take $d=8.122$ kpc [35]111Although the
integration volume might be larger for a more realistic estimation, we limited
the range of the l.o.s. integration to 8.122 kpc as a first step of the study.
This can be justified as a conservative estimation by having confirmed the DM
flux from the Galactic center increase of 40% (80%) if we change the range of
the line-of-sight integration to 10 kpc (38 kpc)..

For the numerical calculation, we use the cosmic ray flux data from the
Galprop codebase [36]. Both the cosmic ray sources and the DM profile would be
denser in the vicinity of the Galactic center in a realistic scenario;
however, following reference [4] we assume the flux of cosmic rays is uniform
throughout the Galaxy. The dependence of CR-DM flux on cosmic ray profiles
near the Galactic Center will be the scope of our future study.

Figure 1: Coordinate system used in the figures showing the numerical results.

  2. 2.

The obtained CR-DM incoming flux in the Galactic rest frame is transformed to
the laboratory frame in the Solar system.

  3. 3.

DM-nucleon scattering caused by the CR-DM in the Solar system is simulated
using the CR-DM flux derived in step 2. Then we obtain the nuclear recoil
energy $E_{R}$ and its scattering angle $\gamma_{\mathrm{G.C.}}$.

We assume three typical density profiles for $\rho_{\chi}(r)$. The first is
the frequently used Navarro–Frenk–White (NFW) profile suggested by N-body
simulations [37],

$\displaystyle\rho_{\mathrm{NFW}}(r)=\frac{\rho_{0}}{\frac{r}{r_{s}}\left(1+\frac{r}{r_{s}}\right)^{2}}.$
(2.2)

The second is the Einasto profile which gives a better fit to the N-body
simulation data than the NFW profile [38, 39]. It is given by

$\displaystyle\rho_{\mathrm{Ein}}(r)=\rho_{0}\exp{\left[-\frac{2}{\alpha}\left\\{\left(\frac{r}{r_{s}}\right)^{\alpha}-1\right\\}\right]}.$
(2.3)

Because both the NFW and Einasto profiles are cuspy in the Galactic center, we
also investigate a cored profile to see the profile dependence clearly. The
third profile is the pseudo-isothermal profile (PIT) [40]

$\displaystyle\rho_{\mathrm{Iso}}(r)=\frac{\rho_{0}}{1+\left(r/r_{s}\right)^{2}}.$
(2.4)

Parameters used in the numerical calculation are shown in Table 1, which are
taken from [35] for the NFW and the Esnasto profiles, and [41] for the PIT
profile. In the table, the local DM density is represented as
$\rho_{\mathrm{DM},\odot}$, and $\rho_{0}$ is normalized to reproduce that
value of $\rho_{\mathrm{DM},\odot}$.

Table 1: Parameters on DM profiles |  $\rho_{\mathrm{DM},\odot}$ [GeV/cm3] |  $\rho_{0}$ [GeV/cm3] |  $r_{s}$ [kpc] | $\alpha$
---|---|---|---|---
NFW | 0.38 | 0.83 | 11.0 | -
Einasto | 0.38 | 0.30 | 9.2 | 0.18
PIT | 0.35 | 3.56 | 2.7 | -

##  3 Numerical Result and discussion

In direct detection experiments with directional sensitivity, the energy and
angle of nuclear recoils are expected to be detected. We will discuss the
results obtained from simulations in this section.

###  3.1 Angular distribution of CR-DM and nuclear recoil

Figure 2 shows the sky maps for the incoming distributions of the CR-DMs of 10
MeV mass expected in the Solar system for the three density profiles assumed.
The left, middle, and right columns correspond to the NFW, the Einasto, and
the PIT profiles, respectively. In the NFW and Einasto profiles, the incoming
flux of CR-DM is highly concentrated in the direction of the Galactic center.
Whereas in the PIT profile, the incoming flux of CR-DM is found over a wide
area centered at the Galactic center. There is almost no difference in the
directional distributions of the CR-DM in the mass range of 1 GeV to 1 MeV.

Next, we move to the directional distributions of nuclear recoils by the CR-
DMs. In Figures 3-5, the angular distributions of the recoiled nuclei are
shown in whole sky map density plots for the following reference cases:

  * •

the targets are $p$, F, and Ag,

  * •

the DM profiles are NFW, Einasto, and PIT,

  * •

the DM masses are 100 MeV, 10 MeV, and 1 MeV.

These figures are intended to illustrate an ideal situation without any
consideration of detector response (energy threshold, energy resolution,
quenching factor and others) and any types of background. The center of the
figures corresponds to the arrival from the direction of the Galactic center,
thus in all simulations, the CR-DM flux from the Galactic center is
significantly higher than the contribution from other directions. In the whole
sky map density plots, there is no noticeable difference between the
simulations using the NFW and Einasto profiles. However, there is a distinct
difference between the simulations using the PIT profile and the other two
profiles. The different appears because, the NFW and Einasto profiles in the
Galactic center are cuspy, but it is cored in the PIT profile. This cuspy vs
cored property is reflected in the fact that events are less dense in the
Galactic center for simulations using the PIT profile than the other
simulations using the NFW and Einasto profiles.

###  3.2 Energy dependence of the angular distribution

We now discuss the distributions of the nuclear recoil energies and scattering
angles to see the dependence of the CR-DM mass, which is not visualized in the
sky maps in subsection 3.1. The differential angular distributions expected
with the $p$, F, and Ag targets are shown in Figure 6, 7, and 8, respectively.
The rate of detection is shown in units of event number per second per
kilogram. NFW, Einasto and PIT profile simulations are shown in the left,
center, and right columns of each figure. Simulations with a CR-DM mass of 100
MeV, 10 MeV, and 1 MeV are shown in the top, center, and bottom rows of each
figure. SI interactions with the DM particles are assumed for the $p$ and Ag
targets, while SD interactions are assumed for the F target. See Appendix A
for the sensitivity with the F target for SI interactions. We assume the cross
section of the DM and the proton interaction is $\sigma_{\chi-p}=10^{-32}$
cm2, which is currently not experimentally constrained for SI or SD
interactions with $m_{\chi}\lesssim 0.1$ GeV [4].

A detector energy threshold of $E_{\mathrm{thr}}=10$ keV is assumed for the
$p$ and F targets, and threshold of $E_{\mathrm{thr}}=100$ keV for the Ag
target is assumed. The event rates expected in the energy range greater than
the detector energy threshold and which elastic scattering occurs are shown in
the figures [42]. The black histograms show the distribution over the whole
energy range of interest, while colored ones show the distributions
corresponding to the recoil energy ranges. We have discarded large momentum
transfer events before making these angular distributions, i.e., events with
energies greater than 10, 3, and 0.6 MeV are discarded for $p$, F, and Ag
targets, respectively. The reason large momentum transfer events are discarded
is that inelastic scatterings are expected to dominate once the kinetic energy
of the DM is large enough to see the internal structure of the nucleus. To
discuss inelastic scatterings between the DM and target nuclei, we would need
to introduce some assumptions depending on the type of DM and its interactions
[43, 44, 45, 11]. In this paper, we prefer to discuss the properties of the
CR-DM in general manners as much as possible, thus we focus only on the
elastic scatterings. We assume that if the de Broglie wavelength of the DM is
longer than that of the target nuclei, it can be treated as a perfect elastic
scattering. For the sake of consistency and considering current detector
setups, only such events are included here, thus we conservatively focus only
on the elastic scatterings.

As can be easily expected, forward scattering events make peaks at
$\cos\gamma_{\mathrm{G.C.}}=1$, and some distributions have peaks in the
center at $\cos\gamma_{\mathrm{G.C.}}=0$. These central peaks are clearly
visible in the distributions for heavy CR-DMs and disappear in ones for light
CR-DMs. The central peaks correspond to the instance where the target and CR-
DM do not change the direction of their motions from their initial ones in the
center-of-mass frame through scattering. This kind of scattering, in which the
target and CR-DM graze each other in the center-of-mass frame, becomes less
visible in the laboratory frame as the CR-DM mass decreases. For light CR-DM,
in the center-of-mass frame, the CR-DM hits a light target and both the CR-DM
and the target are bounced away in opposite directions from its initial state.
As a result, a prominent peak is produced at $\cos{\gamma_{\mathrm{G.C.}}}=1$,
which is the direction expected to be naively scattered from the Galactic
center.

The energy dependence is similar for the simulations of NFW and Einasto
profiles. On the other hand, in the figures there is a difference between the
simulations using the PIT profile and the other two profiles, because of the
different treatment of structures in the center of the galaxy. Now we will
take a closer look at the simulations using the PIT profile. The central peak
near $\cos{\gamma_{\mathrm{G.C.}}}=0$ is moderate and less noticeable than in
the NFW and Einasto profiles. This is because the PIT profile has a core at
the Galactic center. It is clear from Figures 3-5 that the CR-DM signal comes
from a wider range than that in the simulations using the other two profiles,
in which the CR-DM comes from a narrow range near the Galactic center.

###  3.3 Asymmetry and sensitivity

$p$ F
---

Ag
---

Figure 9: Simulated asymmetry for each target and CR-DM mass. Colors indicate
the density profiles. Scattering cross section for DM and proton is assumed to
be $\sigma_{\chi p}=10^{-32}$ cm2, and the exposure is assumed to be NIT 5.0
kg-yr for the target $p$ and Ag, and SF6 1550 kg-yr for the target F.

$p$ F
---

Ag
---

Figure 10: Null-rejection significance of the asymmetry parameter $A$ for each
target. The CR-DM masses are assumed to be 10 MeV and 100 MeV, and are marked
with dotted and solid lines, respectively. The scattering cross section for DM
and proton is assumed to be $\sigma_{\chi p}=10^{-32}$ cm2, which is the same
as the preceding figures.

To quantify the forward-backward asymmetries of the recoil direction
distributions, a parameter $A$ is defined as

$\displaystyle A=\frac{n_{+}-n_{-}}{n_{+}+n_{-}},$ (3.1)

where $n_{+}$ and $n_{-}$ are the number of events for
$\cos\gamma_{\mathrm{G.C.}}>0$ and $\cos\gamma_{\mathrm{G.C.}}<0$,
respectively.

Figure 9 shows the asymmetry parameters and their statistical errors obtained
by the simulation. The colors represent the three density profiles and the
errors are statistical ones. The number of events is normalized as 5.0 kg-yr
exposure of the nano imaging tracker (NIT) for the target $p$ and Ag, and 1550
kg-yr exposure of SF6 gas. The latter is based on the assumption of the
exposure 10k m3 SF6 gas at 20 Torr and six years, which is equivalent to 10
times the exposure envisioned for the Cygnus-1000 project [32]. The weight
ratios of $p$ and Ag in the NIT are set at 1.6% and 44.5%, respectively [46].
For light CR-DMs less than $O(10)-O(100)$ MeV mass, enough statistics clearly
confirm the directionality originated by signals coming from the Galactic
center direction. Because the number of events discarded for the light CR-DM
instance, is smaller than for the heavy CR-DM instance and the central peak
disappears. As already mentioned, events with large recoil energy, which often
occur in heavy CR-DM cases, are discarded because they are not elastic
scattering, thus the number of events available for the analysis is reduced
and the error is significant. The CR-DM generally can have a larger scattering
cross section than WIMPs with masses around O(1 - 100) MeV. The large cross
section often results in attenuation of flux before reaching underground
detectors [4, 5, 6]. In the plots, we assume that the CR-DM descends from the
upper atmosphere and reaches the underground detectors and neglect the
attenuation effect. Due to their large scattering cross section, the component
reaching the detector from the opposite side of the Earth experiences
significant attenuation, thus it is not effective to use this component. Data
selection is feasible for the time intervals when the detector is oriented
towards the central region of the Galaxy, as gas detectors such as the NEWAGE
offer time resolution, and in the NEWSdm, in principle it is possible to
deactivate the nuclear emulsion by temperature adjustment [47]. Thus, if such
an event selection by time is performed, the net exposure would need to be
doubled.

Figure 10 shows the evolution of the null-rejection significance of the
asymmetry parameter $A$ with each target, where the range of the exposure
corresponds up to six years for the same configuration as that in Figure 9.
The CR-DM masses are assumed to be 10 MeV (dotted lines) and 100 MeV (solid
lines). As can be inferred from Figure 9, assuming elastic scattering, all
targets investigated have sufficient potential for the null-rejection of the
asymmetry parameter $A$ within the scope of the envisioned future upgrades.
With enough exposure, the asymmetry can be confirmed with more than 5 $\sigma$
confidence level for 10 MeV (100 MeV) mass CR-DM except for $p$, at most 8
$\sigma$ (7.6 $\sigma$) and 2.6 $\sigma$ (5.2 $\sigma$) for F and Ag,
respectively. For SD interactions using the fluorine target, the NEWAGE
experiment is developing a 1 m3 scale detector filled with SF6 gas.
Furthermore, an update is planned to the CYGNUS-1000 experiment [32], which
has 1000 m3 scale gaseous time projection chamber (TPC). It is planned to
achieve 5 $\sigma$ sensitivity for the 10 MeV CR-DM mass after 2.6 (2.4)
years’ of measurement for Einasto (NFW) profile. For the NEWSdm experiment,
the device production system at the Gran Sasso National Laboratory (LNGS) has
been operated successfully, and it is possible to produce the 10 kg scale
target using the current system. In addition, new scanning systems with a 0.5
kg/year/machine are being developed, and several kg-scale scanning will be
achieved by using the same upgrades on five existing machines. Furthermore, a
high scanning speed machine with a $>$ 5 kg/year/machine has already been
designed that uses a wide field of view with multi-camera imaging and optimal
driving motion system. From the high scanning machine, several 10 kg scale
experiment will be achieved. Low energy threshold tracking will be achieved by
the finer grain nuclear emulsions [46] and super-resolution techniques [48].
Therefore, taking into account the possibility of future upgrades, both NEWAGE
and NEWSdm would have the potential to verify the light CR-DM arriving from
the direction of the Galactic center.

##  4 Conclusion

Up-scattering of DM by cosmic rays is an undeniable key phenomenon that can
allow the search for light DMs, which are not easy to reach by standard
nuclear recoil methods. The CR-DMs are expected to arrive mostly from the
direction of the Galactic center, where the density of the DM is large, and
thus can be examined with direction-sensitive detectors. In this paper, we
study the detection possibilities of CR-DMs with directional detectors and
investigate the expected energy and angular distributions.

We find that almost independent of the DM density profiles, CR-DM tends to
come from the direction of the Galactic center, and the targets in the
detector are scattered towards the forward direction. Assuming elastic
scattering of DM and target, we also investigate the sensitivity of
directional detectors to verify the existence of CR-DM. Our study showed that
the asymmetry of nuclear recoil is expected to hold the characteristic
signature of the CR-DM. It is shown that the feature can be verified by
directional detectors, especially for light DM with mass less than
$O(100)-O(10)$ MeV. With scattering cross section of DM and proton as
$\sigma_{\chi-p}=10^{-32}$ cm2 and $m_{\chi}=10$ MeV ($m_{\chi}=100$ MeV), the
asymmetry can be evaluated at 0.7$\sigma$, 8$\sigma$, and 2.6$\sigma$
(0.2$\sigma$, 7.6$\sigma$, and 5.2$\sigma$) at most for the $p$ target (the SI
interaction), for the F target (the SD interaction), and for the Ag target
(the SI interaction), respectively. For the evaluation, exposures 0.48 kg-yr,
7200 kg-yr, and 13.2 kg-yr are assumed for $p$, F, and Ag targets,
respectively. These large exposures can be achieved with future large
detectors like CYGNUS and NEWSdm.

NFW Einasto PIT
---
Figure 2: Sky map of the CR-DM in the Solar system. The left-handed side, the
center, and the right-handed side correspond to the NFW, the Einasto, and the
PIT profiles, respectively. The CR-DM mass is 10 MeV. The density of the
number of events is written using arbitrary units (A.U.).

NFW, $m_{\chi}=100$ MeV Einasto, $m_{\chi}=100$ MeV PIT, $m_{\chi}=100$ MeV
---

NFW, $m_{\chi}=10$ MeV Einasto, $m_{\chi}=10$ MeV PIT, $m_{\chi}=10$ MeV
---

NFW, $m_{\chi}=1$ MeV Einasto, $m_{\chi}=1$ MeV PIT, $m_{\chi}=1$ MeV
---

Figure 3: Sky map of the nuclear recoil directions in the Solar system. The
left column, the center column, and the right column correspond to the NFW,
the Einasto, and the PIT profiles, respectively. The upper row corresponds to
the CR-DM mass of 100 MeV, the middle row to 10 MeV, and the lower row to 1
MeV. The nuclear target is $p$. We assume SI interactions between the DM
particle and the nuclear target.

NFW, $m_{\chi}=100$ MeV Einasto, $m_{\chi}=100$ MeV PIT, $m_{\chi}=100$ MeV
---

NFW, $m_{\chi}=10$ MeV Einasto, $m_{\chi}=10$ MeV PIT, $m_{\chi}=10$ MeV
---

NFW, $m_{\chi}=1$ MeV Einasto, $m_{\chi}=1$ MeV PIT, $m_{\chi}=1$ MeV
---

Figure 4: Legend is same as Figure 3 except that the target atom in the direct
detection is F.

NFW, $m_{\chi}=100$ MeV Einasto, $m_{\chi}=100$ MeV PIT, $m_{\chi}=100$ MeV
---

NFW, $m_{\chi}=10$ MeV Einasto, $m_{\chi}=10$ MeV PIT, $m_{\chi}=10$ MeV
---

NFW, $m_{\chi}=1$ MeV Einasto, $m_{\chi}=1$ MeV PIT, $m_{\chi}=1$ MeV
---

Figure 5: Legend is same as Figure 3 except that the target atom in the direct
detection is Ag.

NFW, $m_{\chi}=100$ MeV Einasto, $m_{\chi}=100$ MeV PIT, $m_{\chi}=100$ MeV
---

NFW, $m_{\chi}=10$ MeV Einasto, $m_{\chi}=10$ MeV PIT, $m_{\chi}=10$ MeV
---

NFW, $m_{\chi}=1$ MeV Einasto, $m_{\chi}=1$ MeV PIT, $m_{\chi}=1$ MeV
---

Figure 6: Angular distribution of the nuclear recoil in each energy band, in
the Solar system. The left column, the center column, and the right column
correspond to the NFW, the Einasto, and the PIT profiles, respectively. The
upper row corresponds to a CR-DM mass of 100 MeV, the middle row to 10 MeV,
and the lower row to 1 MeV. The nuclear target is $p$. We assume SI
interactions between the DM particle and the nuclear target.

NFW, $m_{\chi}=100$ MeV Einasto, $m_{\chi}=100$ MeV PIT, $m_{\chi}=100$ MeV
---

NFW, $m_{\chi}=10$ MeV Einasto, $m_{\chi}=10$ MeV PIT, $m_{\chi}=10$ MeV
---

NFW, $m_{\chi}=1$ MeV Einasto, $m_{\chi}=1$ MeV PIT, $m_{\chi}=1$ MeV
---

Figure 7: Legend is same as Figure 6 except that the target atom in the direct
detection is F.

NFW, $m_{\chi}=100$ MeV Einasto, $m_{\chi}=100$ MeV PIT, $m_{\chi}=100$ MeV
---

NFW, $m_{\chi}=10$ MeV Einasto, $m_{\chi}=10$ MeV PIT, $m_{\chi}=10$ MeV
---

NFW, $m_{\chi}=1$ MeV Einasto, $m_{\chi}=1$ MeV PIT, $m_{\chi}=1$ MeV
---

Figure 8: Legend is same as Figure 6 except that the target atom in the direct
detection is Ag. As mentioned in subsection 3.2, we assume
$E_{\mathrm{thr}}=100$ keV in the analysis of asymmetry, but data for
$10<E_{R}<100$ keV are also shown here for reference.

## Acknowledgments

It is a pleasure to thank S. Matsumoto, K. Hata, K. Hayashi, C. Cappiello and
H. Oshima for helpful discussions. KIN is also grateful to S. Kime, T. Tomiya
and N. Ikeda for discussion on the early stages of this work. This work was
supported by JSPS KAKENHI Grant Number 19H05806, 26104005, (A) 16H02189, (A)
18H03699, (C) 21K03562, (C) 21K03583, 21K13943, 22H04574 and Wesco Scientific
Promotion Foundation.

## References

  * [1] E. Aprile, et al.  Dark Matter Search Results from a One Ton-Year Exposure of XENON1T.  Phys. Rev. Lett., Vol. 121, No. 11, p. 111302, 2018.
  * [2] J. Aalbers, et al.  First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment.  7 2022.
  * [3] Wen Yin.  Highly-boosted dark matter and cutoff for cosmic-ray neutrinos through neutrino portal.  EPJ Web Conf., Vol. 208, p. 04003, 2019.
  * [4] Torsten Bringmann and Maxim Pospelov.  Novel direct detection constraints on light dark matter.  Phys. Rev. Lett., Vol. 122, No. 17, p. 171801, 2019.
  * [5] Yohei Ema, Filippo Sala, and Ryosuke Sato.  Light Dark Matter at Neutrino Experiments.  Phys. Rev. Lett., Vol. 122, No. 18, p. 181802, 2019.
  * [6] Christopher V. Cappiello and John F. Beacom.  Strong New Limits on Light Dark Matter from Neutrino Experiments.  Phys. Rev. D, Vol. 100, No. 10, p. 103011, 2019.  [Erratum: Phys.Rev.D 104, 069901 (2021)].
  * [7] James B. Dent, Bhaskar Dutta, Jayden L. Newstead, and Ian M. Shoemaker.  Bounds on Cosmic Ray-Boosted Dark Matter in Simplified Models and its Corresponding Neutrino-Floor.  Phys. Rev. D, Vol. 101, No. 11, p. 116007, 2020.
  * [8] Gang Guo, Yue-Lin Sming Tsai, and Meng-Ru Wu.  Probing cosmic-ray accelerated light dark matter with IceCube.  JCAP, Vol. 10, p. 049, 2020.
  * [9] Shao-Feng Ge, Jianglai Liu, Qiang Yuan, and Ning Zhou.  Diurnal Effect of Sub-GeV Dark Matter Boosted by Cosmic Rays.  Phys. Rev. Lett., Vol. 126, No. 9, p. 091804, 2021.
  * [10] Joerg Jaeckel and Wen Yin.  Boosted Neutrinos and Relativistic Dark Particles as Messengers from Reheating.  JCAP, Vol. 02, p. 044, 2021.
  * [11] Gang Guo, Yue-Lin Sming Tsai, Meng-Ru Wu, and Qiang Yuan.  Elastic and Inelastic Scattering of Cosmic-Rays on Sub-GeV Dark Matter.  Phys. Rev. D, Vol. 102, No. 10, p. 103004, 2020.
  * [12] Nicole F. Bell, James B. Dent, Bhaskar Dutta, Sumit Ghosh, Jason Kumar, Jayden L. Newstead, and Ian M. Shoemaker.  Cosmic-ray upscattered inelastic dark matter.  Phys. Rev. D, Vol. 104, p. 076020, 2021.
  * [13] Keiko I. Nagao and Hiroshi Okada.  Modular A4 symmetry and light dark matter with gauged U(1)B$-$L.  Phys. Dark Univ., Vol. 36, p. 101039, 2022.
  * [14] Jie-Cheng Feng, Xian-Wei Kang, Chih-Ting Lu, Yue-Lin Sming Tsai, and Feng-Shou Zhang.  Revising inelastic dark matter direct detection by including the cosmic ray acceleration.  JHEP, Vol. 04, p. 080, 2022.
  * [15] Chen Xia, Yan-Hao Xu, and Yu-Feng Zhou.  Production and attenuation of cosmic-ray boosted dark matter.  JCAP, Vol. 02, No. 02, p. 028, 2022.
  * [16] Xiangyi Cui, et al.  Search for Cosmic-Ray Boosted Sub-GeV Dark Matter at the PandaX-II Experiment.  Phys. Rev. Lett., Vol. 128, No. 17, p. 171801, 2022.
  * [17] R. Xu, et al.  Constraints on sub-GeV dark matter boosted by cosmic rays from the CDEX-10 experiment at the China Jinping Underground Laboratory.  Phys. Rev. D, Vol. 106, No. 5, p. 052008, 2022.
  * [18] Chen Xia, Yan-Hao Xu, and Yu-Feng Zhou.  Azimuthal asymmetry in cosmic-ray boosted dark matter flux.  6 2022.
  * [19] K. Abe, et al.  Search for Cosmic-ray Boosted Sub-GeV Dark Matter using Recoil Protons at Super-Kamiokande.  9 2022.
  * [20] Gilly Elor, Robert McGehee, and Aaron Pierce.  Maximizing Direct Detection with HYPER Dark Matter.  12 2021.
  * [21] Debjyoti Bardhan, Supritha Bhowmick, Diptimoy Ghosh, Atanu Guha, and Divya Sachdeva.  Boosting through the Darkness : Bounds on boosted dark matter from direct detection.  8 2022.
  * [22] James Alvey, Torsten Bringmann, and Helena Kolesova.  No room to hide: implications of cosmic-ray upscattering for GeV-scale dark matter.  9 2022.
  * [23] Tarak Nath Maity and Ranjan Laha.  Cosmic-ray boosted dark matter in Xe-based direct detection experiments.  10 2022.
  * [24] Moqbil S. Alenazi and Paolo Gondolo.  Directional recoil rates for WIMP direct detection.  Phys. Rev. D, Vol. 77, p. 043532, 2008.
  * [25] Ciaran A. J. O’Hare and Anne M. Green.  Directional detection of dark matter streams.  Phys. Rev. D, Vol. 90, No. 12, p. 123511, 2014.
  * [26] Bradley J. Kavanagh.  Discretising the velocity distribution for directional dark matter experiments.  JCAP, Vol. 07, p. 019, 2015.
  * [27] Bradley J. Kavanagh and Ciaran A. J. O’Hare.  Reconstructing the three-dimensional local dark matter velocity distribution.  Phys. Rev. D, Vol. 94, No. 12, p. 123009, 2016.
  * [28] Riccardo Catena, Jan Conrad, Christian Döring, Alfredo Davide Ferella, and Martin B. Krauss.  Dark matter spin determination with directional direct detection experiments.  Phys. Rev. D, Vol. 97, No. 2, p. 023007, 2018.
  * [29] Keiko I. Nagao, Tomonori Ikeda, Ryota Yakabe, Tatsuhiro Naka, and Kentaro Miuchi.  Discrimination of anisotropy in dark matter velocity distribution with directional detectors.  Phys. Dark Univ., Vol. 27, p. 100426, 2020.
  * [30] M. Andriamirado, et al.  Limits on sub-GeV dark matter from the PROSPECT reactor antineutrino experiment.  Phys. Rev. D, Vol. 104, No. 1, p. 012009, 2021.
  * [31] F. Mayet, et al.  A review of the discovery reach of directional Dark Matter detection.  Phys. Rept., Vol. 627, pp. 1–49, 2016.
  * [32] S. E. Vahsen, et al.  CYGNUS: Feasibility of a nuclear recoil observatory with directional sensitivity to dark matter and neutrinos, arXiv:2008.12587.  8 2020.
  * [33] A. Aleksandrov, et al.  NEWS: Nuclear Emulsions for WIMP Search.  4 2016.
  * [34] N. Agafonova, et al.  Discovery potential for directional Dark Matter detection with nuclear emulsions.  Eur. Phys. J. C, Vol. 78, No. 7, p. 578, 2018.
  * [35] P.F. de Salas, K. Malhan, K. Freese, K. Hattori, and M. Valluri.  On the estimation of the local dark matter density using the rotation curve of the milky way.  Journal of Cosmology and Astroparticle Physics, Vol. 2019, No. 10, pp. 037–037, oct 2019.
  * [36] A.E. Vladimirov, S.W. Digel, G. JÃ³hannesson, P.F. Michelson, I.V. Moskalenko, P.L. Nolan, E. Orlando, T.A. Porter, and A.W. Strong.  Galprop webrun: An internet-based service for calculating galactic cosmic ray propagation and associated photon emissions.  Computer Physics Communications, Vol. 182, No. 5, p. 1156â1161, May 2011.
  * [37] Julio F. Navarro, Carlos S. Frenk, and Simon D. M. White.  A Universal density profile from hierarchical clustering.  Astrophys. J., Vol. 490, pp. 493–508, 1997.
  * [38] J. F. Navarro, E. Hayashi, C. Power, A. R. Jenkins, C. S. Frenk, S. D. M. White, V. Springel, J. Stadel, and T. R. Quinn.  The inner structure of $\lambda$cdm haloes - iii. universality and asymptotic slopes.  Monthly Notices of the Royal Astronomical Society, Vol. 349, No. 3, p. 1039â1051, Apr 2004.
  * [39] Alister W. Graham, David Merritt, Ben Moore, JÃ¼rg Diemand, and BalÅ¡a TerziÄ.  Empirical models for dark matter halos. ii. inner profile slopes, dynamical profiles, and rho/sigma3.  The Astronomical Journal, Vol. 132, No. 6, p. 2701â2710, Jan 2006.
  * [40] Raul Jimenez, Licia Verde, and S. Peng Oh.  Dark halo properties from rotation curves.  Mon. Not. Roy. Astron. Soc., Vol. 339, p. 243, 2003.
  * [41] Christopher Wegg, Ortwin Gerhard, and Marie Bieth.  The gravitational force field of the galaxy measured from the kinematics of RR lyrae in gaia.  Monthly Notices of the Royal Astronomical Society, Vol. 485, No. 3, pp. 3296–3316, feb 2019.
  * [42] T. Ikeda, T. Shimada, H. Ishiura, K.D. Nakamura, T. Nakamura, and K. Miuchi.  Development of a negative ion micro tpc detector with sf6 gas for the directional dark matter search.  Journal of Instrumentation, Vol. 15, No. 07, p. P07015, jul 2020\.
  * [43] John R. Ellis, R. A. Flores, and J. D. Lewin.  Rates for Inelastic Nuclear Excitation by Dark Matter Particles.  Phys. Lett. B, Vol. 212, pp. 375–380, 1988.
  * [44] J. Engel and P. Vogel.  Neutralino inelastic scattering with subsequent detection of nuclear gamma-rays.  Phys. Rev. D, Vol. 61, p. 063503, 2000.
  * [45] G. Arcadi, C. Döring, C. Hasterok, and S. Vogl.  Inelastic dark matter nucleus scattering.  JCAP, Vol. 12, p. 053, 2019.
  * [46] Takashi Asada, Tatsuhiro Naka, Ken-ichi Kuwabara, and Masahiro Yoshimoto.  The development of a super-fine-grained nuclear emulsion.  PTEP, Vol. 2017, No. 6, p. 063H01, 2017.
  * [47] M. Kimura, T. Naka, S. Furuya, T. Asada, T. Katsuragawa, M. Yoshimoto, A. Umemoto, S. Machii, H. Ichiki, O. Sato, and Y. Hoshino.  WIMP tracking with cryogenic nuclear emulsion.  Nucl. Inst. Meth. A, Vol. 845, pp. 373–377, 2017.
  * [48] Andrey Alexandrov, Takashi Asada, Di Crescenzo Antonia De Lellis, Giovanni, Valerio Gentile Gentile, Tatsuhiro Naka, Valeri Tioukov, and Atsuhiro Umemoto.  Super-resolution high-speed optical microscopy for fully automated readout of metallic nanoparticles and nanostructures.  Sci. Rep., Vol. 10, No. 18773, 2020.

##  Appendix A The SI and SD cross sections

Fluorine is the target assumed in the directional direct detection through SD
interactions. The sensitivity to SI and SD cross sections varies by a factor
$\eta_{A}=\sigma_{\chi-p}^{\mathrm{SI}}/\sigma_{\chi-p}^{\mathrm{SD}}$. By
evaluating the factor, the sensitivity to SD interactions in Section 3.2 can
be scaled to that to SI interactions. The cross sections of DM and nucleus
scattering for the SI and SD interactions $\sigma_{\chi-N}^{\mathrm{SI,SD}}$
are

$\displaystyle\sigma_{\chi-N}^{\mathrm{SI}}$
$\displaystyle=\sigma^{\mathrm{SI}}_{\chi-p}\frac{\mu_{\chi-N}^{2}}{\mu_{\chi-p}^{2}}A^{2}\mathrm{[cm^{2}]},$
(A.1) $\displaystyle\sigma_{\chi-N}^{\mathrm{SD}}$
$\displaystyle=\sigma^{\mathrm{SD}}_{\chi-p}\frac{\mu_{\chi-N}^{2}}{\mu_{\chi-p}^{2}}\frac{\lambda^{2}J(J+1)}{0.75}\mathrm{[cm^{2}]},$
(A.2)

respectively, where $\sigma^{\mathrm{SD,SI}}_{\chi-p}$ is the DM and proton
scattering cross section, $\mu_{\chi-N(p)}$ is the reduced mass of CR-DM and
nucleus (proton), $A$ is mass number, and $\lambda J(J+1)=0.647$ for F.
Therefore, the ratio of event number required to reach the same cross section
$\sigma_{\chi-p}$, which is also denoted as
$\eta_{A}=n_{\mathrm{SI}}/n_{\mathrm{SD}}$, is estimated as

$\displaystyle\eta_{A}$ $\displaystyle=\frac{0.75A^{2}}{\lambda^{2}J(J+1)}$
$\displaystyle\simeq 418.$ (A.3)